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Preface 


This research brief contains the progress reports of the Research Staff of the 
Center for Modeling of Turbulence and Transition (CMOTT) from May 1991 to May 
1992. It is intended as an annual report to the Institute for Computational Mechanics 
in Propulsion and NASA Lewis Research Center. A separate report entitled, “Work- 
shop on Engineering Turbulence Modeling,” covering some of the 1991 CMOTT 
Summer research activities was released earlier this year. 

The main objective of the CMOTT is to develop, validate and implement the 
turbulence and transition models for practical engineering flows. The flows of interest 
are three-dimensional, incompressible and compressible flows with chemical reaction. 
During this period, the research covers two-equation (e.g., K-e) and algebraic Reynolds- 
stress models, second moment closure models, probability density function (pdf) 
models, Renormalization Group Theory (RNG), Large Eddy Simulation (LES) and 
Direct Numerical Simulation (DNS). Last year was CMOTT’s second year in 
operation. CMOTT now has eleven members from ICOMP, NASA LeRC, and 
Sverdrup Technology Inc., working on various aspects of turbulence and transition 
modeling in collaboration with NASA-Lewis scientists and Case Western Reserve 
University faculty members. The CMOTT members have been continuously and 
actively involved in the international and national turbulence research activities. A 
biweekly CMOTT seminar series has been conducted with speakers invited from within 
and without of the NASA Lewis Research Center, including foreign speakers. The 
current CMOTT roster and its organization are listed in Appendix A. Listed in 
Appendix B are the abstracts and the scientific and technical issues discussed in 
biweekly CMOTT seminars. Appendix C gives a list of references which are the papers 
contributed by CMOTT members in the last two years. 
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Lumley’s Energy Cascade Dissipation Rate 
Model for Boundary-Free Turbulent Shear Flows 


1. Motivation and Objective 

True dissipation occurs mainly at the highest wavenumbers where the eddy sizes 
are comparatively small. These high wave numbers receive their energy through 
the spectral cascade of energy starting with the largest eddies spilling energy into 
the smaller eddies, passing through each wavenumber until it is dissipated at the 
microscopic scale. However, a small percentage of the energy does not spill contin- 
uously through the cascade but is instantly passed to the higher wave numbers 1 . 
Consequently, the smallest eddies receive a certain amount of energy almost imme- 
diately. As the spectral energy cascade continues, the highest wave number needs 
a certain time to receive all the energy which has been transferred from the largest 
eddies. As such, there is a time delay, of the order r, between the generation of 
energy by the largest eddies and the eventual dissipation of this energy. 

For equilibrium turbulence at high Reynolds numbers, there is a wide range where 
energy is neither produced by the large eddies nor dissipated by viscosity, but is 
conserved and passed from wavenumber to higher wavenumber. The rate at which 
energy cascades from one wavenumber to another is proportional to the energy 
contained within that wave number. This rate is constant and has been used in 
the past as a dissipation rate of turbulent kinetic energy. However, this is true 
only in steady, equilibrium turbulence. Most dissipation models contend that the 
production of dissipation is proportional the production of energy and that the 
destruction of dissipation is proportional to the destruction of energy. In essence, 
these models state that the change in the dissipation rate is proportional to the 
change in the kinetic energy. This assumption is obviously incorrect for the case 
where there is no production of turbulent energy, yet energy continues to cascade 
from large to small eddies. If the time lag between the onset on the energy cascade 
to the destruction of energy at the microscale can be modeled, then there will be a 
better representation of the dissipation process. Development of an energy cascade 
time scale equation will be discussed in Section 2.2. 

2. Work Accomplished 
2.1 Mean Flow Equations 

For incompressible flow, the equations for continuity of mass and axial momentum 
are written as 
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where, u^Uj is the turbulent Reynolds stress tensor. Using the eddy viscosity 
concept 2 , the Reynolds stress may be related to the mean strain rate and a turbulent 
eddy viscosity, 

-UiUj = v t 

The turbulent viscosity given in the above equation can be interpreted as a measure 
of the turbulent kinetic energy existing in the flow times a local length scale, 


dUi dUi 

dxj 


dxj 


kfi a • 
3 3 


( 1 ) 


v t 


ck*l 


where c is an arbitrary constant. The definition used for this length scale is the pri- 
mary discriminating factor between turbulence models and determines the number 
of equations which need to be solved. The length scale can be written in terms of 
the turbulent kinetic energy,/:, and its dissipation rate,e, 



e 


Now, the momentum equation can be written as 

DUt _ d (, .dUA 1 OP 

Dt dxj \ U Ut dxj ) p dxi 

where, 


e 

and is a constant. The k and e are determined by solving transport equations 
for the turbulent kinetic energy and the dissipation. These equations and the mod- 
ification to the “standard” dissipation equation will be discussed in Section 2.2. 

2.2 Turbulence Equations 

Taking the Navier-Stokes equations and multiplying by Ui = Ui + Ui and then 
taking the time average yields an equation that contains both the energy equation 
for the mean velocity, Ui, and a mean energy equation for the fluctuating component 
of velocity, Eliminating the energy equation for the mean velocity yields the 
turbulent energy budget in the absence of a pressure gradient field 3 , 


UiUj S{j 2 us ij S{j . (3) 

The last term in equation (3) is the rate at which viscous forces perform deformation 
work to the fluctuating strain rate, defined as 

1 / duj dui \ 

Sl3 2 V dxi dxj ) 


Dk 

~Dt 


_a_ 

dx i 


UiUiUi — 2vUiS 


l^lj 


2 


Lumley’s Energy Cascade Dissipation Rate Model 


This term is the true viscous dissipation rate which drains kinetic energy from the 
system, e = 2i 'SijSij. This dissipation rate is important to the overall turbulence 
structure. The next to the last term in equation (3) serves to transport kinetic 
energy between the mean flow and the turbulence. It is generally referred to as the 
kinetic energy production term. Writing the mean strain, Sij, as 

2 l dx, 9 Xj ) 

the production term for incompressible flow is defined as 


Pk 


U , XI j S{j — IS t 


dUi + dUi\ 

dxi dxj ) 


dUj 

dxj 


( 4 ) 


The last set of terms to be modeled represents the energy redistribution by viscous 
forces and are contained within the parenthesis on the right hand side of equation 
(3). From the definition of the fluctuating strain rate, 


2 vmsij = v 


dk 

dxj 


and using a mean gradient hypothesis, the triple correlation can be rewritten as 


—UiUiUj 


v t dk 

(x k dx j 


where Ck is a constant of the order one. Combining all these modeling assumptions 
yields the high Reynolds number form of the turbulent kinetic equation, 

Dk _ _d_ f / vA dk_\ fdUi dUj\ dUj 
Dt dxj V V + cr k ) dxj / + ^ \ dxj + dxi ) dxj 


An equation of the following form has been widely used for modeling of the 
dissipation, 

De __ d f / v t \ de \ 1 e 2 

Dt dxj \{ U+ a e ) dxj ) € ck 

The production term, P e has been assumed to be proportional to the turbulent 
kinetic energy production term, equation (4). This approximation places a di- 
rect correlation between the production of turbulent energy and the production of 
dissipation. Although this may be true at equilibrium, where by definition the dis- 
sipation equals the turbulent energy production, there are some deficiencies in this 
statement for the more general case. According to Lumley 1 , the modeled dissipation 
rate equation should be written as, 


Wt = ((" + ”) I 1 ) + “ €<s ~ 

Dt dxj \\ o f J dxjj cp cd k 


( 5 ) 
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The production term, P e , now represents the rate at which energy arrives at the dis- 
sipative wavenumbers through the energy cascade starting with low wavenumbers. 
In equation (5), 1/S represents the characteristic time for energy to be transported 
to the dissipative wave numbers. 

In the following equation, the term on the left hand side and the first term on the 
right represent the transport of the inverse time scale by advection and diffusion, 
respectively. The remaining terms are the production from the mean straining 
forces within the large scale mean flow structures and the dissipative losses at the 
microscale, 


PS 

Dt 
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The production and loss terms have been scaled by a characteristic time, T , which 
in this study has been taken to be a weighted function of the time scale for the large 
eddies, r, 

m T 

T = . 

2Cb 


Since most energy enters the cascade at scales comparable to the integral scale, l, 
the characteristic time scale of the large eddies dictates the rate of energy transfer 
and is defined as r = £ for nearly isotropic turbulence. Using the approximation 

for the integral length scale as, l = the time scale becomes 


r 



which is valid even in nonequilibrium flows. Letting Cb be a function of mean and 
turbulent quantities, the inverse time scale equation can be written as 



where 


C B = 1 + 


9 k^/SijSij 


Determination of the modeling constants is discussed in the next section. 


2.3 Model Constants 

In this paper, only high Reynolds number turbulence has been considered. The 
model constants are examined for the grid turbulence case where the mean velocity 
gradients and strain rates vanish. For this case, turbulent intensity decays with the 
following proportionality, 


A 

k o 
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where, t = — 
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so that the dissipation rate must decay as 
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oc 


ar +i) 


Here, n is the decay rate and is typically in the range 1.1 < n < 1.3. Substituting 
these relationships into the inverse time scale equation yields 



raCg 

tc B 


dt. 


After integration, this equation is 


S = at 


>c £ 

C B 


where a is a constant of integration. For there to be an inverse scaling relationship 
between S and t, the exponent in the above equation must reduce to -1 for the grid 
turbulence case. Since the coefficient C B reduces to 1.0 in grid turbulence, then 


c B = ra- 


in this study, n has been set to 1.1 and held constant 4 , which sets the value of c B . 
Now using the fact that, initially, 



the coefficients ca and c B are related as follows 

1 _ 3 c A ra + 1 

cd 2\/2 c D n 

The coefficient c ^ has been evaluated for equilibrium turbulence where, where 
P k = e. Rewriting equation (2) as 

Vtt 

” k 2 


and using equations (1) and (4), this coefficients is 


c M — 


a i Uj 


0.09. 


After the above analysis, the only adjustable parameters are ^ and the a coeffi- 
cients, which should be of the order one. The complete set of optimized coefficients 
is given in Table 1. 


O’ k = 0.9 

n = 1.1 

o e = 1.1 

to 

II 

3 

os = 1.0 

C B = 1 + gk y/ s 'i s a 

= 0.09 

1 3 CA I n+ 1 

CD “ 2v^2 C D ' n 

^ = 1.13 

CD 



Table 1. k-e-S Model Coefficients 
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2.4 Free Shear Flows 

The three turbulent transport equations are solved implicitly using the finite 
difference technique of Spalding 5 for parabolic equation systems. An initial profile 
is imposed upon the first marching plane for the turbulent kinetic energy and the 
dissipation, and <S is initially set to S, = In the mixing layer and jet cases, 
the initial plane was equally divided into a uniform velocity field and still air. For 
the case of the plane wake, this initial profile came from a flat plate calculation 
computed with the same solution technique. Starting with the initial plane, the flow 
field is developed by marching the transport equations downstream and applying 
the following boundary conditions, 

Planar Jet 

Centerline (y = 0) g=0,g=0 

0 _ SoidCkkUy^o+^^iv+vt)^) 

C> V = 0 — c t fc!7 v= o + eAx 

Far Field (y = y max ) fc = 0 ,e = 0,<S = 0 
Round Jet 

Centerline (r = 0) §£ = 0 , = 0 

e S 0 idCbkU r= o + $;((i'+i't)^') 

C>r = (i ~ c b kU r=0 +eAx 

Far Field (r = r TOax ) A: = 0 ,e = 0,<S = 0 
Planar Mixing Layer 

Centerline (y = y m in) 773 —^ — - = 0.015, — = 3.5 x 10 _G 

u y=ymtn kj y _ yrnin 

c — 1 

^y—ymxn — y. 

Far Field (y = y max ) fc = 0 ,e = 0,<S = 0 
Planar Wake 

Centerline (y = 0) S ^=°>^ =0 

0 _ SoldCbkU v ^o + $;((i'+Vt)^') 

<> y=0 — c b fcJ7 v= o + eAx 

Far Field (y = y m ax) k = 0,e = 0,S = 0 

The solutions have been checked after several hundred steps to insure that the 
profiles have become self-similar. 

2.4.1 Planar Jet 

The three equation k-e-S model performed very well for planar jet flow when 
compared to the experimental data of Gutman and Wygnanski 6 , Bradbury 7 and 
Heskestad 8 . Mean velocities profiles given in Figure 1 are well predicted by both the 
fc-e-S model and the standard k-e model, which is reassuring since the coefficients 
for both of these models were optimized for the planar jet case. The spreading rates 
predicted by both models are within the experimental range and are given in Table 
2 . 

There is a slight discrepancy between the k-e-S model and the standard k-e 
model for centerline kinetic energy as can be seen in Figure 2. However, both 
predictions are well within the range of experimental measurements, which show 
a considerable spread near the centerline of the jet. Both predictions are closely 
aligned with the data of Bradbury 7 . 
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Shear stress is also well predicted (see Figure 3) where the standard k-e model 
over-predicts the shear stress at the edge of the jet in comparison to the data of 
Gutman and Wynanski 6 and Bradbury 7 . The k-e - S model predicts a narrower 
profile, also indicated by the smaller spreading rate, which agrees with Bradbury’s 
measurements. Towards the centerline of the jet, both models accurately predict 
the increase in shear stress and the extrema near % = 0.08. 

For the next two figures, Figures 4 and 5, there is no experimental data available 
for comparison. The dissipation reaches a maximum at the location of the peak 
shear stress, Figure 4. From this extrema, the dissipation decreases toward the 
centerline and also eventually trails off to zero at the edge of the jet. Figure 5 
shows the behavior of the time scale for the planar jet. Notice that the minimum 
energy transfer time between the large scale structure to the dissipative microscale 
is in the middle of the jet (where 5 is a maximum since the dimension of S is < -1 ). 

2.4.2 Axisymmetric Jet 

A primary problem with the standard k-e model is its predictive capability for 
a 2-D versus a 3-D jet, the free jet anomaly. In Figures 6 - 10, the standard k-e 
model is unable to predict any of the turbulence quantities correctly, overpredicting 
not only the turbulent kinetic energy and the shear stress but also the mean velocity. 
The standard k-e can be “fixed” using a suggestion by Pope 9 . Writing the standard 
e equation for reference, 


the correct spreading rate can be obtained by changing C\ from its standard value of 
1.45 to 1.6 (see Table 2). Clearly, by modifying a coefficient, the standard k-e model 
is capable of predicting the correct spreading rate; however, this model cannot be 
considered general. The spreading rate predicted by the k-e - S model is very close 
to the experimentally measured rate and has been obtained with no modifications 
to the model. 

Comparisons between the A:-e-S model and the experimental data 10,11 ’ 12 , Figures 
6-10, indicates that the new model accurately predicts the mean velocity profiles 
and the shearing stress across the jet. Less accurately predicted by the k-e - S model 
is the centerline turbulent kinetic energy, Figure 7, where there is considerable 
scatter in the experimental data. 

2.4.3 Planar Mixing Layer 

For the case of a planar mixing layer, neither model accurately predicts the high 
speed side of the layer (see Figures 11 - 15). The experimental data 13 contains 
a much more dispersive edge compared to the computations, Figure 11. One ex- 
planation may be wind tunnel noise. Adding a boundary condition specifying a 
given noise level, 12 percent fluctuating velocity, increases the spreading on the 
high speed side of the computations and improves the comparisons with the data. 
This boundary condition is used with both models. Unfortunately, there is no way 
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to determine what the true noise level was in the experiments or whether its effect 
has been properly accounted for computationally. 

An analytical solution for the mean velocity profile has been described by Schlichting 14 . 
This analytic solution lies between the experimental data and the two computational 
models for the case when a = 9.5 ( a is a “tweeking” constant used to tune the ana- 
lytic solution to match the experimental data. Typical values of a are between 9 13 
and ll 15 .) 

Both the k-e-S and the standard k-e model predict the peak value of turbulent 
kinetic energy quite well (Figure 12). However, this peak is shifted toward the low 
speed side of the mixing layer compared to the experimental peak. Patel 16 indi- 
cated that the data of Wygnanski and Fielder 13 had shifted towards the low speed 
side, which is opposite from the trend noticed in this study. Since the discrepancy 
between the data and the computations lies on the high speed side of the mixing 
layer, the problem could be a result of wind tunnel noise as previously mentioned. 

In the the plot of Reynolds stresses, Figure 13, both models over predict the shear 
stress by nearly 30 percent and consequently over-predict the spreading rate given 
in Table 2. 

The k-e-S model and the standard k-e model predict almost identical levels of 
dissipation (Figure 14). Since a turbulent kinetic energy level is specified at the high 
speed boundary, the dissipation is also specified on this boundary. A constant value 
of v l L - = 6.5 x 10 -6 gave the best fit to the data. Figure 15 shows the inverse time 
scale* distribution through the mixing layer. The minimum energy transfer time, 
where S is a maximum, corresponds to the location of the maximum shear stress. 

At the high speed edge of the mixing layer, the time scale is specified to b e<S=|. 

2.4.4 Planar Wake 

For the case of the planar wake, the mean velocity profile is correctly predicted 
towards the centerline of the wake, but drops off too quickly at the wake edges (see 
Figure 16). Although the mean velocity shows the correct shape in Figure 16, the 
centerline velocity is high compared to the analytical solution described in Reference 

dy ± 

[14] (see Figure 17). As a result, the spreading rates, predicted by 

the turbulence models are 6 percent too low for the k-e-S model and 3 percent too 
low for the k-e model as compared to the analytical solution. In other words, the 
turbulence models are predicting a much more compact wake with a smaller velocity 
defect than is seen in the experimental data or the analytical fit. This trend is also 
seen in the turbulent kinetic energy profiles in Figure 18 where the k-e-S model can 
predict the correct maximum intensity but drops off too quickly towards the edges. 
Also, although the computed peak value for shear stress is correct the width of the 
curve is too narrow, with the location of the extrema lying too near the centerline 
(see Figure 19). 

The dissipation levels computed by the models apparently reduce the turbulent 
intensity of the flow too much near the edges of the wake, which prevents the wake 
from spreading. There is unfortunately no experimental data to compare with the 
dissipation curves in Figure 20. Interestingly, the location of the minimum energy 
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transfer time in Figure 21 is towards the outer edges of the wake where the mean 
flow should exert its greatest influence. 



Experiment 

k-e-S Model 

k-e Model 

Planar Jet 

0.11-0.12 

0.114 

0.114 

Round Jet 

0.085-0.095 

0.095 

0.126, ci = 1.45 
0.0928, ci = 1.6 

Planar Mixing Layer 0.16 

0.186 

0.186 

Planar Wake 

.101 14 

0.095 

0.098 


Table 2. Spreading Rate Comparisons for Free Shear Flow 


3. Future Work 
Future plans include: 

i) Extending this model to wall bounded flow. Specifically, a low Reynolds number 
version of this model will be developed for near wall turbulence. 

ii) Continuing to run this model for more test cases over a variety of flow fields. 
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Figure 1: Mean velocity profile for a turbulent, planar Figure 3: Shear stress profile for a planar jet. U ma 
jet. Umax- centerline velocity. centerline velocity. 



Figure 2: Turbulent kinetic energy profile for a planar Figure 4: Dissipation rate profile for a planar jet. Umax : 
jet. U ma x : centerline velocity. centerline velocity. 
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Figure 5: Time scale profile for a planar jet 



Figure 6: Mean velocity profile for a turbulent, axisym- 
metric jet. Umax : centerline velocity. 


Figure 7: Turbulent kinetic energy profile for an ax- 
isymmetric jet. Umax • centerline velocity. 



Figure 8: Shear stress profile for an axisymmetric jet. 
Umax • centerline velocity. 
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Figure 11: Mean velocity profile for a planar mixing 
Figure 9:. Dissipation rate profile for an ax.symmetnc ^ centerline velocity. — analytical solution 

jet. U m ax • centerline velocity. 



Figure 12: Turbulent kinetic energy profile for a planar 
Figure 10: Time scale profile for an axisymmetric jet layer. Umax • centerline velocity. 
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Figure 13: Shear stress profile for a planar mixing layer. 
Umax ■ centerline velocity. 



Figure 14: Dissipation rate profile for a planar mixing 
layer. U max ■ centerline velocity. 



Figure 15: Time scale profile for a planar mixing layer 



freestream velocity. — analytical solution 
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Figure 17: Axial mean centerline velocity for a planar 
wake behind a flat plate. Uoo'. freestream velocity. 


Figure 19: Shear stress profile for a planar wake be 
hind a flat plate. U cen ur : centerline velocity, t/co 

freestream velocity. 
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Figure 21: Time scale profile for a planar wake behind 
a flat plate 
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PDF Turbulence Modeling and DNS 

A. T. Hsu 

1. Motivations and Objectives 

The problem of time discontinuity (or jump condition) in the coalescence/dispersion 
(C/D) mixing model is addressed in first part (PDF) of this work. A C/D mix- 
ing model continuous in time is introduced. With the continuous mixing model, 
the process of chemical reaction can be fully coupled with mixing. In the case 
of homogeneous turbulence decay, the new model predicts a pdf very close to a 
Gaussian distribution, with finite higher moments also close to that of a Gaussian 
distribution. Results from the continuous mixing model are compared with both 
experimental data and numerical results from conventional C/D models. 

The efFect of Coriolis forces on compressible homogeneous turbulence is studied 
using direct numerical simulation (DNS). The numerical method used in this study 
is an eigth order compact difference scheme. Contrary to the conclusions reached 
by previous DNS studies on incompressible isotropic turbulence, the present results 
show that the Coriolis force increases the dissipation rate of turbulent kinetic energy, 
and that anisotropy develops as the Coriolis force increases. The Taylor-Proudman 
theory does apply since the derivatives in the direction of the rotation axis vanishes 
rapidly. A closer analysis reveals that the dissipation rate of the incompressible 
component of the turbulent kinetic energy indeed decreases with a higher rotation 
rate, consistent with incompressible flow simulations (Bardina 17 ), while the dissipa- 
tion rate of the Compressible part increases; the net gain is positive. Inertial waves 
are observed in the simulation results. 

2. Work Accomplished 

2.1 PDF Turbulence Model for Combustion 

In Collaboration with J.Y. Chen 

Accurate prediction of turbulent reacting flows requires the solution of an evo- 
lution equation for the probability density function (pdf) of the thermo-chemical 
variables using Monte Carlo simulation. Since the pdf equation, like most equations 
describing turbulent motion, is not closed, closure models have to be deyised. For 
the pdf of scalars, the terms in the pdf equation that need modeling are molecular 
mixing and turbulent convection. The present work deals with the modeling of 
molecular mixing. 

Most of the mixing models are based on the coalescence/dispersion (C/D) model 
by Curl 6 This model is known to have deficiencies, and efforts had been made to 
correct these deficiencies, for example, Janicka, et al. 1979 and Pope 1982. The 
most recent efforts have been devoted to the problem of coupling between mixing 
and chemical reaction. Chen and Kollmann 4 proposed a reaction conditioned model 
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that allows correct prediction of combustion in the flame-sheet regime. Norris and 
Pope 10 proposed a new model based on ordered pairing that aimed at the same end. 

All the existing models suffer in one respect, namely, they are discontinuous in 
time: once a pair of particles are chosen to participate in mixing, their properties will 
jump abruptly regardless of the step size of the time integration. This phenomenon, 
clearly non-physical, could cause difficulty in coupling the processes of mixing and 
chemical reaction. In the present work, a new model that is continuous in time is 
proposed. With this new model, the processes of molecular mixing and chemical 
reaction can be fully coupled. 

In the case of homogeneous turbulence decay of a scalar, one expects a Gaussian 
distribution for the pdf, and finite values for the higher moments. Pope 11 pointed 
out the modified Curl model could not produce the correct pdf for this problem, 
and the higher even moments from that model tend to infinity; he suggested an age. 
biased sampling process to overcome these shortcomings. The present continuous 
model, as we will show, predicts a pdf distribution very close to Gaussian for ho- 
mogeneous turbulence decay, and gives finite higher moments with values close to 
that of a Gaussian distribution. 

The continuous mixing model is applied to the study of both non-reacting and 
reacting flows, and the results are compared with earlier calculations by Hsu 7 as 
well as with experimental data. 

2.1.1 Molecular Mixing Model 

The evolution equation of a single point probability density function of scalar ran- 
dom variables i>\ ,...,tp n — representing the species mass fraction and temperature 
— can be written as 

N 

pd t P + pv a d a P + p'Y^d't’i •••> MP} 

1 

= ~d a (p < V«| M X ) = i’k > P) 

N N 

-P (< €i i\M x ) = ipk> p) 

i— 1 j — 1 

where the terms represent the rate of time change, mean convection, chemical re- 
action, turbulent convection, and molecular mixing, respectively; P is the density- 
weighted joint pdf: 

P = pP/P, 

e is the scalar dissipation: 

€ ij — P $ <y(p j 0 ( \(p j i 

(where D is the diffusion coefficient), and < x\y > denotes the mathematical ex- 
pectation of a random function x conditioned upon y. 

The left hand side of the above equation can be evaluated exactly and requires 
no modeling; the right hand side terms contain the conditional expectation of the 
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velocity fluctuation and the conditional expectation of the scalar dissipation, which 
are new unknowns and require modeling. In the present work we concentrate on 
the modeling of the second term, namely, the conditional expectation of the scalar 
dissipation, referred to as molecular mixing in the following. 

The Modified Curl Model 

The simplest and most used mixing model is the modified Curl model, which 
assumes binary interaction between sample fluid particles. As described by Pope 12 
in a Monte Carlo simulation, the continuous pdf is replaced by delta functions 

1 N 
n=l 

where each delta function represents one sample fluid particle of an ensemble of N 
particles. The evolution of P* entails the movement of the particles in the ^-space, 
or the evolution of the individual values of 0 n ’s. 

With the modified Curl model, the change of <j> n due to molecular mixing is 
achieved by the following binary interaction process: divide the flow domain into 
small cells, each containing N sample particles. Given a small time interval 8t and 
a turbulent time scale t, select randomly N mx pairs of particles, 

N mx = 0.5^iV, 

( C = 6.0) and let a pair, say, m and n, mix as follows 

<(>n(t + St) = A<t> m (t) + (1 - A)<t> n (t) ( 1 ) 

<t>m(t + St) = A(j) n (t) + (1 - A)<t>m(t) (2) 

where A = 0.5£, with £ a random variable uniformly distributed on the interval 
[0,1]. The remaining N - 2N mx particles remain unchanged: 

4>n{t + St) = <pn{t) 

This model does not represent the true physical process since the properties of 
the sample particles change discontinuously regardless the size of the time interval 
6t. This deficiency can be best illustrated by rearranging eq. (1) and dividing it by 
6t 

(Mt + it) - M*)) = _ *„(*)) ( 3 ) 

01 ot 

The derivative ^ does not exist because as 6t goes to zero the right hand side of 
the equation becomes infinite since both A and the difference between <pm(t) and 
are finite. This means that there is a sudden jump in the value of the scalar 
quantities, which is typical of a Poisson process, but is non-physical in the present 
case since the flow properties of turbulence are continuous. 
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Continuous Mixing Model 

One can see from the previous section that the modified Curl model relies on the 
parameter N mx to control the extent of mixing. On the individual particle level, 
it assumes complete mixing once the particle is selected as one of the mixing pair, 
without considering the size of 6t. 

In order to achieve continuous mixing, we propose the following model: during a 
time interval 8t, we assume that all the particles within a cell participate in mixing. 
The extent of the mixing is controlled at the individual particle level. That is to 
say, the N particles within a given cell are randomly grouped into N/2 pairs; the 
properties of all the particles change according to eqs. (1) and (2) The extent 
of mixing now has to be controlled at the individual particle level through the 
parameter A, which is redefined as 


where C' = 2.0. With this new definition, eq. (3) can be written, in the limit 
6t -»• 0, 

dt t 

The above equation states that the change of </> n due to mixing is proportional to 
the difference between (j) m and 0„, and inversely proportional to the turbulence 
time scale r. 

The Coupling of Mixing and Reaction 

The processes of mixing and chemical reaction are essentially decoupled when one 
uses the discontinuous C/D models. In contrast, with the above continuous model, 
coupling becomes natural since, for a given particle, mixing and chemical reaction 
can be described with a single equation: 

“XT = — <f>n(t)) + w ni 

dt t 

where w n is the chemical source term. 

Since the continuous mixing model allows full coupling of the reaction and mixing 
processes, the C/D model with reaction zone conditioning by Chen and Kollmann 4 
can be easily implemented in the present model to simulate the fast reaction in the 
flame sheet regime. Here a modified finite difference version of eq. 12 has to be 
used since w n is infinity in case of fast reaction. 

2.1.2 Results and Discussions 

The continuous mixing model described in the previous section has been vali- 
dated using both non-reacting and reacting flow test cases. The results and their 
comparisons with earlier calculations (Hsu 7 ) using the modified Curl model as well 
as with experimental data are presented in this section. 
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Homogeneous Turbulence Decay of a Scalar 

The case of decaying fluctuation of a passive scalar in homogeneous turbulence 
is used to test the continuous mixing model. The initial condition is 

‘AT/ 2 N 

1)+ d(ip + 1) , 

n=l n=AT/2+l 

that is, in the Monte Carlo simulation, half of the particles are ascribed the value 
1, and the other half -1. 

The pdf distribution of the normalized variable (ip— < (p >)/&, where < <p > is 
the mean and a is the standard deviation, in the homogeneous turbulence decay 
problem converges to a single curve after certain time, and the correct distribution 
should be Gaussian. Fig. 1 and 2 are the pdf distributions from the modified Curl 
model and the present model, both compared to a normal distribution. One can 
see that the pdf from the modified Curl model deviates considerable from Gaussian, 
while the result from the present model is fairly close to a Gaussian distribution. 

The evolution history of the rms and fourth and sixth moments of the scalar 
fluctuation are calculated using both the modified Curl model and the present 
continuous model. Fig. 3 shows the results from the modified Curl model. One 
can see that although the rms from that model behaves well, the fourth and sixth 
moments grow quickly out of bound, oscillating at a level several order of magnitudes 
higher than the value of Gaussian distribution. These results are similar to what 
Pope 11 had observed. Fig. 4 shows the results for the same set of quantities from the 
present model. The rms behaves similar to that from the modified Curl model. The 
fourth and sixth moments, on the other hand, are quite different from those of the 
previous model; they rise smoothly to the value predicted by Gaussian distribution. 
Although the values do not seem to converge, they remain finite, and are of the 
same order of magnitude as that of the Gaussian distribution. 

The above results clearly demonstrated the advantage of the present model over 
that of the modified Curl model. Pope 11 had devised an age biased scheme that 
achieved the same end, which required an additional variable, namely the age of 
the particles, and two extra adjustable parameters. In contrast, the present model 
needs no extra work or parameters. 

Heated Turbulent Jet 

Extensive experimental results for a heated turbulent plane jet have been reported 
by many authors (Bashir, et al. 2 , Browne et al. 3 , Uberoi and Singh 13 , Jenkins 9 , An- 
tonia, et al. 1 ). The turbulent jet has a slightly higher temperature than the ambient. 
Measurements of both the the mean temperature and the rms of temperature fluc- 
tuations were given. We compared the solutions for the heated turbulent jet from 
the new model with experimental data as well as with previous solutions (Hsu 7 ) 
obtained using the modified Curl model. 

In the present study, a combined CFD — Monte Carlo algorithm is used. The 
mean flow field is obtained by solving the Navier-Stokes equation and a two-equation 
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turbulence model using a finite difference scheme. The temperature is treated as a 
conserved scalar and simulated by the pdf equation. 

Fig. 5 shows the comparisons of the the mean temperature distribution from 
the pdf Monte Carlo simulations and experimental data from various authors. The 
figure shows that both mixing models predict the mean temperature distribution 
accurately. 

The results for standard variation, or rms, of the temperature distribution are 
given in Fig. 6. Although the two solutions do not show significant difference, the 
new model seems to agree slightly better with the experimental data. The skewness 
and the flatness, i.e., the third and fourth moments of the temperature fluctuation, 
are given in Figs. 7 and 8; the comparisons show that in the present case of a 
turbulent jet the statistical behavior of the new model is similar to that of the 
modified Curl model. 

Hydrogen-Fluorine Diffusion Flame 

The continuous model has been applied to the H 2 —F 2 diffusion flame studied 
earlier by Hsu 7 . The flow conditions are set according to an experiment performed 
by Hermanson and Dimotakis 6 . The flame consists two streams. The upper stream 
contains 96% of N 2 and 4% of F 2 , the flow velocity is U\ = 22 m/s; the lower stream 
contains 96% of N 2 and 4% of H 2 , with velocity U 2 = 8.8 m/s. The estimated 
Damkohler number ranged from 25 to 130 (Hermanson and Dimotakis 6 ), and a fast 
chemistry model is deemed appropriate in the calculation. Again a modified version 
of eq. 12 is used to accommodate the fast chemistry. 

Fig. 9 shows the temperature rise due to combustion. In the figure, 6t is the 
shear layer thickness determined by 1% of the temperature rise, AT is the actual 
temperature rise due to combustion, ( the two streams have the same temperatures 
initially,) and A T a df is the adiabatic flame temperature assuming complete reaction. 
Details on the flow conditions can be found in Hermanson and Dimotakis 6 . The 
agreement between numerical predictions and experimental data is fairly good, and 
a comparison of the results from the continuous model and that from the modified 
Curl model shows that both performed well for this case. 

Combination with Reaction Zone Conditioning 

Chen and Kollmann 4 developed a mixing model based on reaction zone condi- 
tioning, aimed at the coupling of reaction and mixing. We have shown in Section 2.3 
that with the present model, the processes of reaction and mixing can be fully cou- 
pled; therefore it is only natural to apply the reaction zone conditioning suggested 
by Chen and Kollmann here. 

The H 2 — F 2 diffusion flame problem is reformulated such that the chemical reac- 
tion is confined to a very narrow zone near stoichiometry. By applying reaction zone 
conditioning to the continuous mixing model, we were able to produce a scatter plot 
of the temperature vs. mixture fraction in which all the points reached the equi- 
librium temperature. This result is shown in Fig. 10. The mixture fraction here is 
defined as the molar concentration of fuel divided by the total molar concentration, 
and stoichiometry is located at / = 0.5. 
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2.1.3 Conclusions 

A turbulence mixing model that is continuous in time has been introduced. The 
deficiency of non-physical jump condition in the mixing process is removed in the 
new model. It has been shown that the new model is superior to the existing 
modified Curl model (Janicka, et al. 8 ) in that it can predict a Gaussian distribution 
and finite higher moments in the case of homogeneous turbulence decay; it has 
accomplished what the age biased sampling scheme (Pope 11 ) is designed for, without 
the extra parameters required by that scheme. The numerical results from the 
present model compare well with experimental data. 

2.2 DNS of Homogeneous Compressible Turbulence in a Rotating frame 

Rotation is an important factor in many flow phenomena in nature and in engi- 
neering. Problems that are strongly affected by rotation include flows in turboma- 
chinery, large scale motions in the atmosphere and oceans, and galactic motions. 
Turbulence is important in all the above examples. In order to model turbulence 
in a rotating frame, we need a better understanding on the effects of rotation on 
turbulence. 

Several experiments had been carried out by various researchers 14-16 . However, 
the conclusions are inconsistent. Some results show that rotation increases the 
turbulence dissipation rate, others suggest the opposite. 

Numerical studies of turbulence in a rotating frame are few. In fact the only 
two are incompressible flow studies, by Bardina, et al. 1 ' and by Speziale et al. 18 , 
using the same computer program. Studies other than the present one dealing with 
compressible flows in a rotating frame are not known to the authors. 

Bardina et al . 17 performed both large eddy and full simulations of incompress- 
ible isotropic turbulence. Their results show that, in the case of incompressible 
turbulence, the dissipation rate decreases with increasing rotation rate, and that 
anisotropy does not develop as a result of the Coriolis force. Speziale et al. con- 
firmed these results using a smaller Rossby number, i.e. a faster rotation. 

In general, compressible turbulence in a rotating frame is not homogenous or 
isotropic because of the existence of centrifugal force and density fluctuations. How- 
ever, when the rotation rate is low, or when the solution domain is very close to the 
rotation axis, the centrifugal force can be neglected and the flow is approximately 
homogeneous. In the present study, our goal is to identify the effects of the Cori- 
olis force, as opposed to the effects of the centrifugal force or the combined effects 
of the centrifugal and Coriolis forces. Therefore, the centrifugal force is dropped 
from the governing equations, regardless of whether or not it is negligible. A future 
study will concentrate on the non-homogenous effects of the centrifugal force and 
the combined effects of the two forces. 

The numerical results presented in the following show that although the Coriolis 
force does not appear explicitly in the turbulence kinetic energy equation, its effect 
can be significant. Anisotropy does develop as a result of the Coriolis force. The 
dissipation rate of the compressible and incompressible modes of the turbulent ki- 
netic energy behave differently under Coriolis forces, and the combined effect is a 
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increased dissipation rate. 

2.2.1 Method 

The governing equations for compressible flows in a rotating frame are 
P,t + ( P u j),j = 0 

+ (p u j u i),j = P,t + T i j,j ~ P^ijk^klm^j^l x m ~ ^PUjk^jUk 
(p e ),t + (pUje)j = —pujj -f UijTij — (kTj)j 
p = pRT 

The last two terms in the second equation represent, respectively, the effects of 
centrifugal force and the Coriolis force. The centrifugal force is a nonuniform body 
force which would induce nonuniform pressure and density distributions - , its effect 
on turbulence is to destroy homogeneity. The Coriolis force, on the other hand, is 
a uniform body force, provided that the velocity field is uniform. In the case of 
incompressible flows, the centrifugal force can be included into a modified pressure 
term and thus does not appear explicitly; however, this can not be done in the case 
of compressible flows. 

In order to identify the homogeneous effect of the Coriolis force alone, and to 
compare results with incompressible flow simulations, we simply take away the 
centrifugal force in our study. This, of course, greatly simplifies the problem: With 
a homogeneous initial condition, the flow field remains homogeneous under the 
Coriolis force. The simplified problem is a hypothetical one and can not normally 
be found in nature except in very restricted situations such as the ones mentioned 
in the introduction. 

The governing equations are solved using a compact difference scheme. Lele 19 has 
shown that such schemes have spectral accuracy and are well suited for DNS. The 
spatial derivatives are approximated by 8th order finite differences, and the time 
integration is performed using a three-stage Runge-Kutta time marching scheme. 

The initial condition for the rotating flow simulation is obtained from a simulation 
of isotropic decaying turbulence. The initial condition for the decaying turbulence 
is generated using a given 3D spectrum of the form 


After the decaying turbulence shows the correct skewness, the Coriolis force is 
imposed onto the flow field. 

2.2.2 Results and Discussion 

As a primary study, three cases have been calculated, with rotation in the z- 
direction: The nondimensional rotation rates are D = 0, 40, 80, which correspond 
to turbulent Rossby numbers of oo, 0.0075, and 0.00375. These Rossby numbers 
are much smaller than those considered in References 17 and 18. 
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Figures 11 and 12 show the development of a ID two-point correlation and a ID 
spectrum in the turbulent flow field, where square symbols are initial conditions and 
lines represent later times. Figure 13 shows a comparison of the energy spectrum 
obtained using two different grids, one double the size of the other. The agreement 
affirms sufficient resolution. 

Figures 14-19 show the time history of various statistics of the three cases, where 
time is nondimensionalized using the initial eddy turn over time. 

Figure 14 is the ratio of the turbulent kinetic energy, where ql is the initial 
turbulent kinetic energy and q 2 the value at later times. It is immediately clear 
that the Coriolis force causes a faster decay of the turbulence, a conclusion that is 
contrary to what was observed in the incompressible flow simulations by previous 
researchers 17,18 . Nevertheless, as we shall see in the following, the present results 
are in fact consistent with those of Refs 4 and 5. 

To prove that the present results are consistent with the previous incompressible 
simulations, we look at the ratio of the means of the divergence squared and vorticity 
squared, which are proportional to, respectively, the compressible and incompress- 
ible turbulent kinetic energy. Figures 15 and 16 show that while the dissipation 
rate of the compressible part of the turbulent kinetic energy increases with larger 
Coriolis forces, the dissipation rate of the incompressible part decreases. The com- 
bined effect, as we have seen from Figure 14, is an increase in the total dissipation 
rate. Since the simulations of References 4 and 5 are for incompressible flows, it is 
not surprising that only a decrease in dissipation rate was observed. 

In both References 17 and 18, no anisotropy was observed. Bardina et al. 17 
speculated that a reorganization to two-dimensional flows could occur at higher 
rotation rates, while Speziale et al. suggested that Taylor-Proudman reorganization 
would not occur in a rapidly rotating isotropic turbulence. The present results show 
definite signs of anisotropy as soon as the Coriolis force is turned on. The anisotropy 
can be clearly observed from the time history of the normal components of the 
Reynolds stress, given in Figures 17, 18, and 19 for the three rotation rates. With 
the Coriolis force, the x- and y-components of the normal Reynolds stress decrease 
much faster than the z-component. In fact, it seems that the increased dissipation 
is primarily in the x- and y-directions, and the z-component of the turbulent kinetic 
energy is not affected much. 

The Taylor-Proudman theory stipulates that under strong rotation, with the 
rotation axis parallel to, say, the z axis, then for any flow variable u , we have 
(, du/dz ) = 0. Figures 20 and 21 are the contours of the x-component of velocity 
plotted on three planes. Without the Coriolis force (Figure 20), the flow struc- 
ture remains isotropic. With the Coriolis force, columns soon appear in the flow 
field, suggesting that a Taylor-Proudman reorganization has occurred. Using larger 
Rossby numbers such as those used in References 17 and 18, the results show no 
appreciable anisotropy, hence there appears to be no contradiction. 

Finally, inertial waves can be observed from the time history of various statistics 
presented in Figures 14-19. The frequency is approximately Q/7T, the theoretical 
intrinsic frequency for rotating flows. 
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3. Future Work 

3.1 PDF 

The pdf model is being implemented into an existing compressible flow solver: 

the RPLUS2D code. Work on implementing pdf to RPLUS3D will follow. 

3.2 DNS 

To understand the effect of compressibility, incompressible homogeneous turbu- 
lence in a rotating frame will be studied. 
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Figure 1. Asymptotic pdf distribution for a scalar in homo- 
geneous turbulence. — modified Curl model; G an ss i an. 
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Figure 2. Asymptotic pdf distribution for a scaIat in homo- 
geneous turbulence. — present model; Gaussian. 



Figure 3. Evolution of moments from the modified Curl 

modeL — standard deviation, 0.01 x fourth central 

moment, - - 0.0001 x sixth central moment, o 0.01 x fourth 
moment for Gaussian distribution, □ 0.0001 x sixth moment 
for Gaussian distribution. 



Figure 4. Evolution of moments from the present modeL 

standard deviation, 0.1 x fourth central moment, - - 

0.01 x sixth central moment, o 0.1 X fourth moment for 
Gaussian distribution, □ 0.01 x sixth moment for Gaussian 
distribution. 
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V/V.sr. 

Figure 5. Mean temperature in heated plane jet. — contin- 
uous mixing model, - - - modified Curl model, A Browne et 
al. (1984), □ Bashir and Uberoi (1975), o Uberoi and Singh 
(1975), o Jenkins and Goldschmidt (1973). 



Figure 8. Flatness of temperature variance in heated plane 
jet. — continuous mixing model, modified Curl model. 
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Figure 6. RMS of temperature variance in heated plane jet. 

— continuous mixing model, modified Curl model, o 

Antonia et al. (1983), o Bashir and Uberoi (1975), A Uberoi 
and Singh (1975). 



Figure 7. Skewness of temperature variance in heated plane 
jet. — continuous mixing model, - - - modified Curl model. 



Figure 9. Temperature rise in an H 7 -F 7 difusion flame. 
— continuous mixing model, - - - modified Curl model, o 
Hermanson and Dimotakis (1989). 
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Figure 10. Joint pdf between temperature and mixture frac- 
tion on the centerline of the Hr~F 7 diffusion flame; reaction 
is restricted to a narrow zone at stoichiometry; results are 
obtained by applying reaction zone conditioning (Chen and 
Kollmann, 1990) to the continuous mixing model. 
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Fig. 14 Turbulent kinetic energy Fig. 17 Reynolds stresses , 0 =0.. Re=20 




Fig. 1 5 Compresseble energey ratio. 



Fig. 1 6 Reynolds stresses, Q =40.. Re-20 
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Fig. 16 Incompresseble energey ratio. 



Fig. 19 Reynolds stresses, 0 =00., Re=20 
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Figure 20 . Flow structure without Coriolis Force, 
u-velocity. 
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Figure 21 . Flow structure with the Coriolis force, 
u-velocity. 
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Modeling of Turbulence and Transition 

Tsan-Hsing Shih 

1. Motivation and Objective 

The first objective is to evaluate current two-equation and second order closure 
turbulence models using available direct numerical simulations and experiments, 
and to identify the models which represent the state of the art in turbulence mod- 
eling. 

The second objective is to study the near-wall behavior of turbulence, and to 
develop reliable models for an engineering calculation of turbulence and transition. 
The third objective is to develop a two-scale model for compressible turbulence. 

2. Work Accomplished 

2.1 Evaluation of two-equation models (N.J. Lang and T.-H. Shih 1 ) 

Twelve two-equation models including k — e, k — r, k — to and q — u> models, 
have been evaluated using a common flow solver code (GENMIX) for wall bounded 
turbulent flows. For each model, calculations were carried out for two-dimensional, 
fully developed channel and flat plate boundary layer flows. These flows are ap- 
pealing for model testing because they are simple and self-similar, yet demonstrate 
important features of wall bounded turbulent shear flows. In addition, we can 
compare these calculations with Direct Numerical Simulations (DNS). A list of the 
models which were tested are shown in the table below: 


Ch 

Chien 2 1982 

k - e 

Sh 

Shih 3 1991 

k — e 

LB 

Lam and Bremhorst 4 1981 

k - e 

NH 

Nagano and Hishida 5 1987 

k — e 

NT 

Nagano and Tagawa 6 1990 

k — e 

LS 

Launder and Sharma 7 1974 

k - 6 

JL 

Jones and Launder 8 1973 

k — e 

MK 

Myong and Kasagi 9 1988 

k — e 

YS 

Yang and Shih 10 1991 

k — e 

WI1 

Wilcox 11 1984 

k — oj 

WI2 

Wilcox 12 1991 

k — u) 

SAA 

Speziale, Abid and Anderson 13 1990 

k - t 

Co 

Coakley 14 1983 

q — u> 


Two dimensional channel flow calculations were made at Re T = 180 and Re T = 
395. These cases were compared with DNS data of Kim et al 15 . Calculations for 
the two-dimensional flat plate boundary layer flow at Ree = 1410 were compared 
with DNS data of Spalart 16 . Some flat plate boundary layer comparisons were made 
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between the experimental data of Klebanoff 17 at Re$ = 7700 and solutions of the 
various models. The detailed results are reported in NASA TM 105237, 1991. 

An important criterion for two-equation model comparisons is not only how well 
the model predicts mean velocity and shear stress, but also the turbulent kinetic en- 
ergy and dissipation (or specific dissipation) rate. These predictions should provide 
appropriate turbulent velocity and length scales so that the model can be applied 
to more complex flows for which a simple mixing length model often fails. Some 
researchers maintain that it is not critical whether or not the turbulent kinetic en- 
ergy and the turbulent length scale are predicted with great accuracy. However, 
one may imagine that a two-equation model making unreasonable turbulent veloc- 
ity and length scale predictions would be very questionable when applied to more 
general flows. A model which accurately predicts the shear stress and mean veloc- 
ity does not imply that it has correctly modeled the turbulent kinetic energy and 
length scale. In fact, if the turbulent kinetic energy is incorrect, then the length 
scale must also be incorrect to compensate for the error in the turbulent kinetic 
energy. For this case, two wrongs are making a right. This warrants some caution 
when computing flows for other geometries. 

Comments on two-equation models: 

It is clear that the JL, LS, WI1 and WI2 models underpredict the near wall 
turbulent kinetic energy compared to the other models. 

The standard k — e model has been proven to provide good results in the high 
Reynolds number range. It is therefore attractive for a near wall k - e turbulence 
model to approach the standard k - e model away from the wall. The LB, LS and 
YS models are the only k - e models in this study which possess this characteristic. 

Because boundary layer and channel flows are self-similar, the solutions should 
be independent of the initial conditions. However, some of the models (SAA, Co, 
and LB) have difficulty when the initial conditions contain large gradients. The Co 
Model is particularly dependent on the initial conditions. Even slight perturbations 
of the initial conditions will yield noticeably different solutions with this model. 

JL, LS, WI1 and WI2 are the only models which do not contain y+. Damping 
functions which contain y + are not desirable because y + is erroneous near flow 
separations and not well defined near complicated geometries. Unfortunately, these 
are the same models which poorly predict the near wall turbulent quantities. 

The Wilcox models (WI1 and WI2) suffer from a numerically awkward boundary 
condition for uj at the wall: 

- esp 38 y + - 0 

We cannot define u» at y + = 0. We have tried two ways to approximate w as y + 
approaches 0. First, we chose a large number for u> wa u and, second, we used an 
asymptotic u wa u = ^>. Test cases showed that the solution does not converge 
as u>t ua ii — i - oo when y + —* 0 for either model. In addition, both Wilcox models 
underpredict the turbulent kinetic energy peak value for both boundary layer and 
channel flows. 
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model the new unknowns Tij , n^- and At the level of the second order closure, 
these new unknowns are usually modeled with algebraic equations in terms of the 
second moments and the mean quantities (with the exception of the trace e kk — 2e, 
which is modeled with a transport equation). 

In a general turbulent shear flow with moderate inhomogeneity, the turbulent 
diffusion terms in the second moment equations are usually smaller than the other 
terms. However, the pressure-strain rate and dissipation rate tensors are always 
among the leading terms. Therefore, the performance of modeled equations largely 
depends on the models of pressure-strain rate tensor and dissipation rate tensor. 

In this study, we only concentrate on the models of the pressure-strain rate tensor 
and the dissipation rate tensor for the velocity field: n^- — which are modeled 
as 

n a - ea = n<'> + n|f - 

where 11^, 11^ are respectively referred to as the rapid term and the return-to- 
isotropy term. 

Models for the rapid term 

Launder, Reece and Rodi (LRR ): 18 

% = 0.2 Sij + ~-^ 6 (b ik S jk + bj k Si k - hijbuSu) 

2 q 2 22 3 ( 2 . 2 . 1 ) 

jq 7C*2 

H (bi k Qj k -H bj k ^l{ k ^ 

where C-i = 0.4, and 

_mu] 1 

°%j — ^ 

Sij = 2 + Uj,i), 

^ ij — 2^ l ,j ~ Uj,i) 

This model is linear in the Reynolds-stress. It contains only one model constant C 2 . 
This model satisfies the conventional model constraints 34 . It is the most general 
form at the level of linear dependence on the Reynolds-stress. However, as Lumley 24 
pointed out, this model may violate realizability as the turbulence approaches a 
two-component state. 

Speziale, Sarkar and Gatski (SSG ): 19 


n” 1 0.8 - c" „ cjP 

2? 4 ,J 2 ? 1 

H {pik$ jk “1“ bjkSik T^bijbklS kl) 

C 

+ —£-(bi k Clj k + bj k Cli k ) 
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where, 

C* z * = CJ(M«) 1/2 , P = —vTujUij 
C* = 1.8 C 3 * = 1.3, C 4 = 1.25, C 5 = 0.4 

This model is quasi-linear in the Reynolds-stress, because the coefficients in the 
first two terms are not constant, but depend on the invariant of the Reynolds- 
stress tensor and the production P. This model contains four model constants 
(Cf, C3*, C 4 ,C 5 ), therefore one may imagine that it will be difficult to correctly 
calibrate them. In addition, this model does not satisfy the normalization constraint 
which is one of the basic model constraints 34 . If we impose this constraint, then 
the four coefficients will reduce to only one, and this model will reduce to the LRR 
model. Finally, like the LRR model, the SSG model may also violate realizability. 
Fu. Launder and Tselepidakis (FLT): 20 

n ( . 1} 2 

—^Lr = 0.2 Sij + 0.3 (b ik Sjk + bjkSik - - 6 ijb k iS k i ) 

2 q 2 6 

13 

+ (bik^jk "H bjkdik') 

+ 0.2 (b^Sji + b 2 iSu - 2bkjbuSki - 3bijb ki S k i ) (2.2.3) 

+ 0.2 (bh^lji + b 2 iClu) 

+ f’[4&mi(5ifcfl.jfc d" bjkSik) 

+ 126 m i5 n j'(6 Tn fcfInfc "t" b n k^lmk)] 

where r = 0.7, b 2 j = bi k b k j. 

This model is cubic in the Reynolds-stress. The final form selected contains 
one model constant. This model only satisfies a part of the realizability condition, 
corresponding to a two-component state of turbulence. However, when a scalar 
field is involved, this model cannot satisfy Schwarz’ inequality between velocity and 
temperature. This part of realizability is sometimes called joint realizability. 

Shih and Lumlev (SL)i 21 

n ( . 1} 2 

— =- = 0.25jj + 3Ct;j(bikSjk 4" bjkSik ~~ X b{jb k iS k i ) 

2 q 2 6 

+ \(2 - 7 ar 0 )(b ik n jk + b jk Q ik ) (2.2.4) 

+ 0 . 2 ( 6 ^ Sji + b 2 j L Su — 2bkjbuSki ~ 3 bijbkiS^i ) 

+ 0.2 {blfilji + b 2 ji^ln) 

where, 

a 5 = —(1 + 0.8F 1/2 ), F = 1 + 9bijb jk b k i - 
10 ^ 

This model is quasi-quadratic in the Reynolds-stress, because the model coefficient 
a 5 is a function of the invariants of the Reynolds-stress tensor. We emphasize that 
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this model is obtained from a more general form of the expression than the FLT, 
and satisfies both the two component condition and Schwarz’ inequality between 
the velocity and scalar fields. In addition, the final form is simpler than the model 
of FLT. 

Shih and Mansour (SM): 22 


n (1) 2 

— = 0.2 Sij + 3cts(bikSjk + bjkSik — -zbijbkiSki) 

2 q 2 & 

+ — (2 — 7a$)(bikttjk + bjk^ik) (2.2.5) 

o 

+ 0 . 2 (bfiSji + b 2 j t Su — 2 bkjbuSki ~ SbijbkiSki) 

+ 0.2 + b 2 iQ.ii) 

where, £*5 = ^{1 4- 3.5[1 - (1 - F) 1 / 4 ]}. 

This model has the same form as the SL model. It was derived in a different way 
and contains a different model coefficient a 5 which was calibrated from one of the 
DNS data (Rogers 30 ). This model, like the SL model, fully satisfies realizability 
conditions. 

( 2 ) 

Models for the return-to-isotropy term 11^' 

Rotta: 23 

= -eCbij (2.2.6) 

where, C = 3.0. 

This model is linear in the Reynolds stress, and contains one model constant. It 
was widely used and adopted in the LRR model. We notice that this model will not 
allow the turbulence toreach a two-component state, because when any turbulent 
component reduces to q 2 / 9, the model Eq.(2.2.6) will force it to grow. 

Lumley : 24 

ng> = -el/36., + 7 ( 6 ?, + 2//6„/3)] (2.2.7) 

where, 7 = 0 and 


0 = 2 + - exp(— 7.77/\/Re){72/v / Re + 80.1 ln[l + 62.4(-/J + 2.3///)]} 



This model is quasi-linear in the Reynolds stress, because 7 is set to zero, and (3 
is a function of the invariants of Reynolds stress tensor. This model is simple, and 
satisfies realizability. 

Sarkar and Speziale (SS): 25 

l4 2) = -e[C,6;, - 3(C, - 2)(6?, - !&««)] (2.2.8) 
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where C\ = 3.4. 

This is a quadratic model in the Reynolds-stress tensor. It satisfies what is 
called the weak realizability condition. Like the Rotta model Eq.(2.2.8), this model 
will not produce unphysical results. However, it will not allow the turbulence to 
approach a two-component state, which could occur in some situations, for example, 
in near-wall turbulence. 

Haworth and Pope (HP): 26 

n‘f = -e{Cyb,i - C^bij + b% - bUbij + V3)]} (2.2,9) 


where C\ = 8.3, C 2 = 14.8. 

Eq.(2.2.9) is a slow part of the Haworth and Pope’s model for the situations 
with no mean velocity gradient. This model, like the SS model, will not produce 
unphysical results; however, it will also not allow the turbulence to approach a 
two-component state. 

Choi and Lumley (CL): 27 
If III > 0, 

= -etfbij + 7(6^ + 2JWy/3)] (2.2.10.1) 

where, 


(3 = 2 + 


* j?l/2 


7 = 


l + G x 1 

p*F 1 / 2 G 


1 + G x 2 i 

t = (/u/2) 1 / 3 , 


n = (-/// 3) 1 ' 2 


X = ~, G = -x 4 + 0.8x 6 

V 

p' = exp[-9.29/Be‘/ 2 ]{(^ + ^) - [296 - 16.2(* + l) 4 ]//} 
Re=f-~, 11= - bijbij/2 , III = b i:i b jk b k i/3 


If III < 0, 


ng* = Eq.( 28) 


( 2 . 2 . 10 . 2 ) 


The model coefficients in Eq.(2.2.10.1) were obtained using their comprehen- 
sive measurements of turbulence relaxing from axisymmetric expansion. Both 
Eq. (2. 2. 10.1) and Eq.(2.2.10.2) satisfy realizability; however, Eq. (2. 2. 10.1) is valid 
only for III > 0, because £ is not defined when III < 0. 

Craft and Launder (CfcL): 28 


n = -C,s [26,, + 4 Ci(b% - blJijl 3)] - 2 Ebij (2.2.11) 
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where, 

Ci =3.1(A 2 A) 1 ' 2 , C[= 1.2 

9 

A 2 — 4 bijbji) j4>3 — 6 jbjkbkii A 1 — -^ 3 ) 

This model is tensorially quadratic in the Reynolds stress, and satisfies realizability. 
Yamamoto and Arakawa (YA ): 29 

+ Q2 (b% - HUfii/3)] (2.2.12) 

where, 

<*i=2 +p F [q ( b\ k ) r + \bl k \ s sign(bl k )} 

& 2=3 (<*i — 2) 

p = —12, q = —0.65, r = 0.4, s = 0.45 

F = l-^L+96j t 

The YA model tried to fit situations with both positive and negative b 3 kk . It also 
meets the requirement of realizability. 

Shih and Mansour (SfeM ): 22 

ng> = -e{(2.0 + CfFtyn + 7 [b\j + (1/3 + 2 II)bij + |//$y]} (2.2.13) 


where, 

C f = (1/9) exp(-7.77//Re){72//Re + 80.1 ln[l + 62.4(-// + 2.3///)]} 

7 = 70(1-/™), Ae= — 

F = 1+911 + 3 III 
II = — 2 bijbij, III — -bijbjkbpi 
Z = 17/20, V = 1/20, 70 = -2 

This model matches the data of Comte-Bellot and Corrsin 21 and meets the require- 
ment that there will be no return to isotropy in the zero Reynolds number limit. 
This model also satisfies realizability. 

Concluding remarks 

We notice that the Reynolds number in all these simulations is low and therefore 
may not represent real turbulence in nature. However, the model terms concerned 
here are mainly pressure related correlations. The fluctuating pressure is not di- 
rectly related to the viscosity, hence the pressure related correlation terms may 
not be directly affected by the Reynolds number, especially the rapid term. The 
return-to-isotropy term which includes the deviatoric part of the dissipation rate 
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tensor may have some dependence on the Reynolds number. According to the above 
consideration, we think that direct comparisons with low-Reynolds DNS data are le- 
gitimate, although we should keep in mind the possible low-Reynolds number effect 
of the DNS data. 

We have directly compared five rapid models with fifteen DNS flows: four of 
Rogers et al.’s 30 shear flows, eleven of Lee et al.’s 31 irrotational strain flows (axisym- 
metric contraction, axisymmetric expansion and plane strain). Comparing the per- 
formance of the LRR and SSG models, which are tensorially linear in the Reynolds 
stress, we conclude that the SSG model gives very little improvement over the LRR 
model. In fact in many cases, it is worse than the LRR model. The reason is not 
very clear. However, we notice that the SSG model does not satisfy the normal- 
ization condition which may be a cause for its poor behavior. If we impose this 
constraint on the SSG model, then it will exactly reduce to the LRR model. In fact 
it can be shown that the most general form of the rapid model, which is tensorially 
linear in the Reynolds stress, is the LRR model. Therefore, in general, the treat- 
ment used in the SSG model would hardly give any improvement over the LRR 
model. A natural way to improve the model is to use a more general nonlinear form 
and more general model constraints. A typical example is the SL 21 model. It starts 
with the most general form, using full realizability constraints together with the 
other conventional constraints 34 . The result is a well behaved model. Indeed, from 
the direct comparisons with the DNS data, the SL 21 model and its variation form 
of SM 22 model give the best performance in most of the cases. As to the FLT 20 
model, it is also a nonlinear model. It starts with a tensorially cubic dependance on 
the Reynolds stress with constant coefficients (in general, these coefficients should 
not be restricted to constants). In addition, the two-component conditions of tur- 
bulence have been imposed. However, the FLT model ignores Schwarz’ inequality. 
Its final form contains two undetermined model constants, but one of them is set 
to zero. The performance of the FLT model, from the direct comparisons with the 
DNS data, is better overall than the linear models, but does not compare with the 
performance of the SL and SM models. So from these direct comparisons of the 
rapid models, we conclude that the SL 21 model and its variation form SM 22 are 
clearly the best. Having said this we notice that, as Reynolds 33 pointed out, any of 
these rapid models will not show any effect of rotation on the invariants (II, III) of 
the anisotropy tensor bij. This is clearly a theoretical deficiency of the current rapid 
models. A further investigation is needed to find out how serious this deficiency 
will be in practice. 

We have directly compared eight return-to- isotropy models with thirty four DNS 
flows: four shear flows and thirty relaxation flows from axisymmetric contraction, 
axisymmetric expansion and plane strain. As was discussed earlier, all the return- 
to-isotropy models are variations of Eq.(2.2.7) derived by Lumley 24 . Therefore 
the differences in the models are due to the different choices of the model coeffi- 
cients. Two linear models are due to Rotta 23 and Lumley 24 (which is quasi-linear 
in bij). Lumley’s model satisfies realizability, matches the data of Comte-Bellot 
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and Corrsin 32 and the limit of the final period of decaying turbulence. It per- 
forms perfectly when III < 0. It also compares well with the DNS data in which 
III > 0. Rotta’s model does not compare with the performance of Lumley’s model. 
In fact, the nonlinear models of SS, YA, HP and C&L also do not compare with the 
performance of Lumley’s model. Apparently the model coefficients chosen in these 
models are not appropriate. The CL 27 model is designed for flows with III > 0 and 
is based on their experiments on relaxing turbulence. It does work better than Lum- 
ley’s model when III > 0. Finally, the S&M 22 model is a nonlinear model; it works 
just like Lumley’s model when III < 0. When III > 0, it shows an improvement 
over Lumley’s model according to the DNS data. So from these direct comparisons 
of the return-to-isotropy models, we conclude that the combination of Lumley’s 
model and Choi’s model, that is the CL 27 model, will give the best performance. 
The S&M 22 model seems as good as the CL model according to these comparisons. 
Having said this, we notice that the existing return-to-isotropy models do not follow 
the relaxation flows from expansion and plane strain very well. Therefore there is 
still a need to further investigate and improve the return-to-isotropy models. 

For detailed comparison in each flow see the reference 34 . 

2.3 Near-wall behavior of turbulence 

The near-wall behavior of turbulence is re-examined in a way different from that 
proposed by Hanjalic and Launder 35 and followers 36 ' 37 ’ 38 ’ 3 . It is shown that at a 
certain distance from the wall, all energetic large eddies will reduce to Kolmogorov 
eddies (the smallest eddies in turbulence). All the important wall parameters, such 
as friction velocity, viscous length scale, and mean strain rate at the wall, are 
characterized by Kolmogorov microscales. According to this Kolmogorov behavior 
of near-wall turbulence, the turbulence quantities, such as turbulent kinetic energy, 
dissipation rate, etc. at the location where the large eddies become “ Kolmogorov ” 
eddies, can be estimated by using both direct numerical simulation (DNS) data 
and asymptotic analysis of near-wall turbulence. This information will provide 
useful boundary conditions for the turbulent transport equations. As an example, 
the concept is incorporated in the standard k-e model which is then applied to 
channel and boundary layer flows. Using appropriate boundary conditions (based 
on Kolmogorov behavior of near-wall turbulence), there is no need for any wall- 
modification to the k-e equations (including model constants). Results compare 
very well with the DNS and experimental data. 

Here we only list the results from this study, for the detail see NASA TM 105663. 

Model equation and boundary condition 

The K-e equations for incompressible flows can be in general modeled as: 

~ = \(v + —)*•,(],( + VrUijiUij + Uj,i) - £ ( 2 . 3 . 1 ) 

L)l (Tk 

^ = [(<- + -kd,i + iAVij + Uj, i)-| 

Dt c r e .it 

— + vvrUijkUijk ( 2 . 3 . 2 ) 
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C x = 1.44, C 2 = 1.92 

h = 1, h = 1 - 0.22 exp (-^-), R t = ^ ( 2 - 3 - 3 ) 

(j fc = 1, a E = 1.3 

These equations are used only for the flow field outside of the turbulence limit point 
y v , where is non-zero. Therefore, Eq.(2.3.2) will not have singularity problems 
and will not need any near-wall modifications like other K-e models do. 2,8 
Eddy viscosity: 

vr = CJ*. (2.3.4) 

where, 

= 0.09 

fn = [1 — exp(ai.Rfc + a^Rf. + asi?^)] 2 

a, = -1.5 * 10 -4 a s = -1.0 « KT* o 5 = -5.0 * 10" 10 < 2 ' 3 ' 5 ) 

„ - Rl,2 y 

V 

Boundary conditions: at y v = 6v/u T , 

£„ = 0.251 -— (2.3.6) 

1 V 

K v = 0.250« 2 (2.3.7) 

In practical applications, R er and R eoo are large numbers, hence y n lL (L is the 
length scale of a flow field) is usually very small. Therefore, as an approximation 
we may let y^/L = 0, but e v and K v must be given by Eqs.(2.3.6) and (2.3.7) 
respectively. These equations have been applied to the calculations of channel and 
boundary layer flows. 

Comparison of models 

To compare the present model with the DNS data and other models (e.g. Jones 
and Launder 8 , and Chien 2 ), we have made calculations on two channel flows 15 ’ 40 
and two boundary layer flows 16,17 . In the present model, all the model constants 
are the same as used in the standard K-e model 39 . Therefore the present model 
will also be suitable for flows away from the wall. The other two models used here 
for comparison do not have this property. Results are shown in figures 1-4. In 
figure 1 and figure 2, three models are compared with two DNS data for channel 
flows: one with R eT = 180, the other with R eT = 395. The profiles of mean velocity, 
Reynolds stress, turbulent kinetic energy and its dissipation rate are plotted in these 
figures. The present model is significantly better than the other two models. Figure 
3 shows the similar comparison for a turbulent boundary with R e0 = 1410. The 
agreement between the present model and DNS data is excellent. Figure 4 shows 
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the results compared with Klebanoff 17 and other boundary layer experiments. The 
skin friction from DNS data 4 * * * * * * * * * * * 16 is also shown in this figure. The results of present 
model are more consistent with the DNS data than the experiments. 

It is also worthwhile to emphasize that the present model equations with the 
standard model coefficients have the simplest form among all two-equation models. 
Hence, we expect that they will have less numerical stiffness in complex turbulent 
flows. 

2.4 Modeling of transition (Z. Yang and T.-H. Shih) 

A model of intermittancy based on the shape factor is added to a two-equation 
k-e model for prediction of boundary layer transition with a free stream turbulence. 
The detailed model and calculations are give by Yang in this year’s Research Briefs. 

2.5 Modeling of compressible turbulence (W.W. Liou and T.H. Shih) 

A two-scale model is proposed based on Hanjalic-Launder’s multiple-scale concept 
for compressible turbulence, in which a distinct scale created by the compressibility 
is modeled separately by considering the effects of pressure-dilatation and dissipa- 
tion dilatation on large-scale energy transfer rate. The detailed model is given by 
Liou in this year’s Research Briefs. 

2.6 Direct numerical simulation of compressible flows (A. Hsu 
and T.-H. Shih) 

In order to have a better understanding of the effect of compressibility on tur- 
bulence, especially the effect of the formation of eddy-shocklets on turbulence, a 
direct numerical simulation of compressible homogeneous shear turbulent flows is 
been performing. The data of all turbulence statistics are very useful for turbulence 
modeling. The detailed simulation is given by Hsu in this year’s Research Briefs. 

3. Future Plans 

Development of second order closure models: pay special attention to the effects 
of inhomogeneity, non-local property, frame-rotation, compressibility, near-wall be- 
havior in the buffer and log-layers. 
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Figure 1. Comparison of models with DNS of 2-D channel with R er = 
180 . 


46 


Tsan-Hsing Shih 



7 + 7 + 



Figure 2. Comparison of models with DNS of 2-D channel with R er = 
395 . 
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Figure 3. Comparison of models with DNS of 2-D boundary layer flow 
with R e o = 1410. 
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Figure 4. Comparison of models with the experiments of 2-D boundary 
layer flows (Klebanofft 13 ! and others). 
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Turbulence Modeling and Experiments 

Aamir Shabbir 
1. Motivation and Objective 

The best way of verifying turbulence models is to do a direct comparison between 
the various terms and their models 1 ,2,3 . The success of this approach depends upon 
the availability of the data for the exact correlations (both experimental and DNS). 
The other approach involves numerically solving the differential equations and then 
comparing the results with the data. The results of such a computation will depend 
upon the accuracy of all the modeled terms and constants. Because of this it is 
sometimes difficult to find the cause of a poor performance by a model. However, 
such a calculation is still meaningful in other ways as it shows how a complete 
Reynolds stress model performs. 

In this study thirteen homogeneous flows are numerically computed using the 
second order closure models. We concentrate only on those models which use a 
linear (or quasi-linear) model for the rapid term. This, therefore, includes the 
Launder , Reece and Rodi 4 (LRR) model; the isotropization of production 4 (IP) 
model; and the Speziale, Sarkar and Gatski 5 (SSG) model. The purpose of this 
study is to find out which of the three models performs better and what are their 
weaknesses, if any. 

The other work reported here deals with the experimental balnces of the second 
moment equations for a buoyant plume. Despite the tremendous amount of activ- 
ity toward the second order closure modeling of turbulence, very little experimental 
information is available about the budgets of the second moment equations. Part 
of the problem stems from our inability to measure the pressure correlations. How- 
ever, if everything else appearing in these equations is known from the experiment, 
pressure correlations can be obtained as the closing terms. This is the closest we 
can come to in obtaining these terms from experiment, and despite the measure- 
ment errors which might be present in such balances, the resulting information will 
be extremely useful for the turbulence modelers. The purpose of this part of the 
work reported here was to provide such balances of the Reynolds stress and heat 
flux equations for the buoyant plume. 

2.0.0 Work Accomplished 

2.1.0 Comparison of Second Order Models in Homogeneous Flows 

Before presenting the results a note about the LRR model constants used in the 
present study is in order. These constants have evolved to slightly different values 
than those orginally recommended by LRR 4 . The value of the Rotta constant C\ 
(in the return to isotropy term) used in the present study is 3.6 (note that due to a 
different definition of b{j used here the value of C\ differs by a factor of two). The 
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rapid term constant C' 2 was assigned a value of 0.4 in the original LER model. In 
the present study the value used for this constant is 0.55 which is slightly higher 
than the value of 0.5 recommended by Morris 6 . It was found out that the value of 
0.55 led to improvement in the performance of LRR model in all the flows tested 
here. (The improvements were slight for the irrotationally strained flows but 

Figure 1 compares the development of Reynolds stresses computed using these 
three models in a flow through axisymmetric contraction with the DNS data 7 . 

Here we show a typical case of S' = 100.00 ( Sk 0 /e 0 = 55.73, case AXM). All 
the models deviate from the DNS data. However, LRR model gives slightly better 
results than the SSG model with IP model performing the worst. 

Figures 2 and 3 show a similar comparison for flow through axisymmetric expan- 
sion for two different strain rates. For the smaller strain rate flow (S = 0.717, SkJ e Q = 
.408, case EXO) SSG model reproduces thejx 2 development quite well while both IP 
and LRR models underpredict it. For the v 2 component all the models give similar 
results. Therefore, for this low strain rate flow SSG model is better than the other 
two rhodels. For the flow with higher strain rate (5 = 7.17 ,Sk 0 /e 0 = 4.08, case 
EXP)' the LRR model is in excellent agreement with the DNS data for both the 
components while both IP and SSG models show overprediction So for this flow 
LRR model works the best. 

Now we show comparisons for the distortion of turbulence by plane strain for 
four cases of differing strain rates. We start from the lower strain rate case. Figure 
4 compares the evolution of the three non-zero Reynolds stress components for the 
flow with strain rate S = 2.6 ( Sk 0 /e 0 = 2.309, case PXC), For ^component all 
the models underpredict the DNS data. LRR model is slightly better than the SSG 
model. IP model is the worst of the three^For v 2 component IP model works the 
best. LRR model slightly underpredicts v 2 while SSG overpredicts it. The third 
component w 2 is overpredicted by all the models with LRR model being better than 
the other two. Figure 5 shows the similar comparisons for the highest strain rate 
case (5 = 25.0, S 0 k/e 0 = 22.227, case PXE). All the three models underpredict the 
■^component. IP modal is the worst of the three models. LRR model_gives slightly 
better result than t(i|SSG model for this stress component. For v 2 component 
LRR model is the blst and SSG model is the worst of the three. For the w 2 
component all the three models overpredict the DNS data with LRR mdoel being 
closest to the data. 1 From the above four plane strain flow comparisons, we note 
that the performance of all the three models deteriorates as the strain rate increases. 
However, on the overall LRR model works better than the other two models. 

Figure 6 shows the same compariscm with the homogeneous shear flow experiment 8 
(5 = 46.8, Sk 0 /e 0 = 6.46). For the ^component LRRmodel gives the best result 
whereas SSG and IP models overpredict it. For the v 2 component also the LRR 
works the best. SSG model slightly overpredicts the data where as IP model if off 
by a larger margin. For the w 2 component both SSG and IP models reproduce 
the data very well whereas LRR model overpredicts the data. For the shear stress 
component LRR performs reasonably whereas SSG model slightly overpredicts the 
data and IP model is off the data by a higher margin. So, for this experiment, LRR 
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model has better overall performance than the other two models. 

Last, we discuss the evolution of q 2 for the case of rotating homogeneous shear 
flow. Since no experimental or DNS data is available for this flow the comparisons 
will be made (for two cases) with the LES 9 . Bardina 10 pointed out that in this 
case we should be careful in interpreting the comparisons for anything more than 
the trends shown by the LES. In all the cases shown here the initial conditions 
corresponded to isotropic turbulence with eo/Sk 0 = 0.296. Figure 7 shows the 
comparisons for the three cases of different Rosby numbers (= Cl/S). For Cl/S = .25 
we note that all three models significantly underpredict the LES results for q 2 ; SSG 
being closest to the LES data and the LRR being the furthest. Qualitatively all 
the three models reproduce the LES trends. For the case of Cl/S — 0.50 SSG is in 
excellent agreement with the _LES results. Both IP and LRR give identical results 
and give a smaller value of q 2 than the LES. It should be pointed out that SSG 
model constants were partially calibrated against this flow. For the third case of 
Cl/S = 1.0, all the three models give identical results. Since no LES results are 
available for this case the only purpose of showing the results is to see how the 
three models compare with each other. 

2.1.2 Conclusions 

Results were shown from numerical computation of various homogeneous turbu- 
lent flows using three different turbulence models. All of these models use a linear 
(or quasi-linear) model for the rapid part of the pressure strain model. Based on 
their overall performance it is found that LRR model works better than both SSG 
and IP models. For the irrotationally flows the differences between the models and 
DNS data increased with the strain rate with LRR model performing better than 
the other two models. For the simple homogeneous shear flow LRR model better 
than the SSG model (for the DNS both performed equally good but for the exper- 
iment LRR worked better). For the homogeneous shear flows both SSG and LRR 
model showed trends similar to those shown by LES with SSG performing better 
than the LRR model. It is worth noting that SSG model has seven empirical con- 
stants as compared to two in LRR model and on the overall it still does not perform 
better than LRR model. Part of the reason for this may be due to the fact that 
the SSG model does not satisfy the normalization constraint where as LRR model 
does. (Normalization is an exact property of the pressure strain correlation; see 
references 4 and 11 for details.) As has been pointed out by Shih and Lumley 3 , for 
a model of the rapid pressure strain part which is linear in the anisotropic tensor 
and satisfies all of its exact properties, LRR is the most general model. 


2.2.0 Experimental Balances for the Second Moments for a Buoyant Plume 

2.2.1 Heat Flux Budgets 

The transport equation for the vertical (streamwise) heat flux can be written as 


U 


dwt 

dr 


+ W 


dwt 

dz 


i d 


d 


dT — ndT 


- — (ruwt) - —(wwt) - uw— w 2 - 

r dr dz dr dz 


53 


A amir Shabbir 


dr dz p oz 


(i) 


Note that the molecular term is written in local cartesian coordinates. The balance 
of this equation is shown in figure 8. Advection term is the smallest ir^this balance 
and, therefore, contributes least to the transport of the heat flux wt. It is clear 
that in the central core of the flow.(r/z < 0.04), the production of this heat flux is 
maintained by the mean buoyancy gradients and the turbulent buoyancy force i.e. 
the source of energy is the gravitational field. The shear production is relatively 
small in this region. Then there is an intermediate region where the production 
from mean velocity and gravitational field are of the same order. However, for 
r/z > 0.1 (which approximately corresponds to the plume half width), most of the 
production is maintained by the mean velocity and buoyancy gradients and the 
turbulent buoyancy production is only a small fraction of these two. The closing 
term in the heat flux balances is labelled as and represents the sum of the pressure 
correlation and the molecular destruction terms i.e. 


n, = - ( i/ + r) 9tti 8t 


dxi 


dxj dxj 


( 2 ) 


The molecular term in (2) is thought to get weaker with increasing Reynolds and 
Peclet numbers, eventually approaching a value of zero in the limit of local (small 
scale) isotropy. This term was not measured and, therefore, its magnitude relative 
to others can not be established. However, in turbulence modeling, it is customary 
to combine this term with the pressure correlation term® and, therefore, from that 
point of view not knowing each term separately does not reduce the usefulness of 
these budgets. Notice that the shape of this term is very similar to the shape of the 
heat flux wt and its magnitude remains large throughout the flow field. 

The equation for the radial heat flux is 


ug+wg 
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( 3 ) 


The balance of this equation is shown in figure 9. Again, we note that the advection 
term is quite small as compared to the other dominant terms in the equation. Unlike 
the wt heat flux balance, the shear production is extremely small here. This is 
because the gradients of mean radial velocity are much smaller than the gradients in 
the mean vertical (streamwise) velocity. There is no turbulent buoyancy production 
in this equation and all the production is due to the mean buoyancy gradients. We 
note that the term representing sum of the pressure correlation and the molecular 
destruction makes up a substantial part of the budget and its shape is similar to 
the radial heat flux. We also note that this budget can not be divided into any 
subregions, where some phenomenon are more dominant than others, because the 
relative magnitude of each of the terms in equation (3) remains the same across the 
flow field. 
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2.2.2 Reynolds Stress Budgets 

The transport equation for the Reynolds stress, within Bussinesq approximation, 
is 


U k ('U’i'U'j ') ,fc — (,‘U’i‘Uk Uj,k T UjUk Ui,k') 

+ piujt + PjUj + i (uiPj + Ujpj) - 


(4) 


where the viscous diffusion term has been neglected since it will be small as com- 
pared to the turbulent diffusion. 

For reasons of convenience, turbulence modelers do not model the pressure cor- 
relation term in the form as it appears in the above equation but re-write it in a 
different form by separating it into a deviatoric and a non-deviatoric part. Two 
ways of doing this have been suggested in the literature and we will look at both of 
these before deciding which one to use in the present study. The traditional way of 
writing this term is 4 


+ Ujp,i) =^p(uij + Ujj) - ^(puiSjk + puj8 ik ),k 


(5) 


where the first term on the right hand side is the deviatoric part. The second term 
is the so called pressure diffusion term. Lumley 12 (1975) has instead suggested the 
following separation 

-(tt*Pj + UjP,i) = - + UjPl) ~ ~ (6) 


where the term in the square brackets is the deviatoric part and the last term on 
the right hand side is the pressure diffusion term. Regardless of which separation is 
employed a correction or model has to be used f or th e correlation pu k . The model 
used here is due to Lumley 8 is given by puk = —q 2 u k /5. This study indicates that 
the use of this model with (5) produces so much pressure diffusion that it negates 
the velocity diffusion (i.e. due to u TufUk). On this basis it was concluded to use the 
separation given by (6) in the present study. (For further details see Shabbir 1 3). 

Therefore, using (6) the equation for the Reynolds stress can be re-written as 


r 2 1 
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where e = eu. Note that anisotropic part of the dissipation part has been combined 
with the pressure correlation term 11 . The term in the curly parenthesis has a zero 
trace and will be denoted by in the rest of the paper. It is this term whose models 
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have been proposed. Note that the above equation is exact since no approximations 
have been used so far. Now we introduce the model for the pressure diffusion term, 
as given above, and with this approximation the above equation becomes 

Uk ( UiUj),k ~ — + Y5 k ~ i^i^k Uj, k + UjU k Ui, k ) + PiUjt + fijUit 

+ $ij — ~ ( 8 ) 

Note that due to the model for the pressure diffusion term this is no longer an exact 
equation and ~ has been used to emphasize this fact. It is this equation which will 
be balanced out with the experimental data and the term will be obtained as 
the closing term. It should be reminded that in addition to the measurement errors, 
any uncertainty in the approximation of the pressure diffusion will also be lumped 
into $ij. 

The equation for the streamwise Reynolds stress w 2 is given by 
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The balance of this equation is shown in figure 10. Advection is the smallest of 
all the terms. Diffusion term is a gain near the center of the plume and a loss 
in the rest of the flow. Also, its magnitude near the center is comparable to the 
other dominant terms in the balance. We note that the buoyancy production is 
comparable to the production due to mean velocity gradients near the plume center 
but over the rest of the flow field the shear production is much larger than the 
buoyancy production. It is also interesting to note that the buoyancy production 
and dissipation rate approximately balance each other. The closing term in this 
balance is and represents the sum of the pressure correlation term and the 
anisotropic part of the dissipation. This term is a loss for the u 2 budget and we 
note that beyond rfz = 0.08 this term and shear production approximately balance 
each other. 

The equation for the radial component u 2 is given by 
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and its balance is shown in figure 11. Obviously the advection of u 2 has the same 
form as the advection of w 2 . The production due to velocity gradients is a loss 
near the plume center and is a gain after about rfz = 0.04. This is because the 
radial gradient of the radial mean velocity is positive near the plume center. The 
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mechanical production term is not large. The diffusion term is a loss over most of 
the flow field and becomes a gain toward the outer edge of the flow field. The sum 
of the pressure correlation term and the anisotropic part of the dissipation rate is 
obtained as the closing term in the budget and represents a gain for u 2 . We further 
note that beyond r / z = 0.08 it approximately balances the dissipation rate. 

Finally we look at the budget for the shear stress uw as shown in figure 12. Its 
equation is given by 


Both advection and the turbulent buoyancy production are of very small magnitude 
and over most of the flow field these approximately balance each other. Neglecting 
these two terms would not cause any significant change in the shear stress balance. 
We note that the diffusion term is not negligible in this budget. The term $ rz is 
essentially balanced by the difference between the shear production and diffusion 
processes. The shape of $ rz is obviously similar to that of the shear stress and its 
peak approximately corresponds to the peak in the shear production. 

3. Future Plans 

3.1 Turbulence Modeling (with T.-H. Shih) 

(a) . Compare the performance of the various non-linear second order models in 
different homogeneous flows in order to find out their strengths and weakneses. This 
will be an extension of the work presented in section 1 of this brief. 

(b) . To develop and test models for turbulent diffusion terms in the Reynolds 
stress equations using Lumley’s theory of third moments 11 . 

3.2 DNS of Bypass Transition (with T.-H. Shih and G. Karniadakis). 

The bypass transition is an important engineering problem due to its relevence 
to turbomachinery environment and, therefore, there is a considerable interest both 
at LeRC and at CMOTT to study this phenomenon. We are interested in carrying 
out the DNS for this problem both in order to provide a data base for the modeling 
efforts of bypass transition at CMOTT and to study its physics. For the former we 
are interesting in finding out what kind of global parameters, if any, are linked to 
the transition process. For the later we are interested in finding out, for instance, 
what is the effect of anisotropy in the free stream turbulence velocity and length 
scale on the transition process. 

These simulations will be designed after the experiments of Sohn and Reshotko 14 
who studied the bypass transition over a flat plate with differing free stream turbu- 
lence intensities. The results of DNS will be compared with these experiments. 
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3.2.1 Numerical Scheme. 

Currently we are exploring the possibility of using a spectral element code for its 
suitability for doing such a DNS. We are inclined to use a spectral element method 
because of its higher accuracy and its ease of local grid refinement. 

The numerical scheme used in the code involves fractional time discretization 
which results in three sets of semi discrete equations. In the first step advection term 
is handled explicitly using a third order Adams-Bashforth scheme. In the second 
step Poisson equation for pressure is solved implicitly and continuity is satisfied. In 
the third fractional step the diffusion terms are accounted implicitly by a second 
order Crank-Nicholson method. 

In order to carry out the spatial discretization the flow domain is first decom- 
posed into macro elements. Each of these macro elements uses a local cartesian 
mesh by employing Gauss-Labatto collocation points. Then within each macro ele- 
ment the flow variables are represented as tensor product of Chebychev polynomial. 
These representations of the flow variables are then substituted into the governing 
equations and discrete equations are obtained by applying the weighted residual 
technique. 

3.2.2 Test Cases to be run 

Several test cases will have to be run in order to validate the code before a full 
DNS can be carried out. First of these is to solve the laminar boundary layer flow 
over a flat plate in order to insure that the numerical method gives the Blasius 
solution. This will also help us explore the various boundary conditions which can 
be used at the top boundary and at the outflow and latter can be used for the 
mean flow during the DNS. After this has been successfully accomplished the most 
unstable mode disturbances based on the linear stability theory will be intorduced. 
This will allow comparing their growth rates (in the linear region) with the solutions 
from the linear stability theory. The third case will be that of suction and blowing 
through the flat plate and the resutls will be compared with those obtained by 
previous workers. 
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Figure 1. Comparison of the models for flow through axisymmetric contraction with the 
DNS data of of Lee and Reynolds (1985). 
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Figure 2. Comparison of the models for flow through axisymmetric expansion with the 
DNS data of of Lee and Reynolds (1985), case LXO. 
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Figure 3 .. Comparison of the models for flow through axisymmetric expansion with the 
DNS data of of Lee and Reynolds (1985), case EXP. 





Figure 4 . Comparison of the models for distortion by plane strain with the DNS data of 
Lee and Reynolds (1985), case PXC, S-2.6. 





Figure 5 . Comparison of the models for distortion by plane strain with the DNS data of 
Lee and Reynolds (1985), case PXE, S=25.0 
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Figure 6 . Comparison of the models for homogeneous shear flow with the data Tavoularis 
and Corrsin (1981). 
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Figure; 7.. Comparison of the models for the rotating homogeneous shear flow with the 
large eddy simulation of Bardina ct al. (1983). 
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Figure 8 . Budget for wl equation. f igure 9 Budget for uw equation. 




Figure 10 Budget for ut equation. Figure 11. Budget for uw equation. 



Figure 12. Budget for u > 2 equation. 
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Modeling of Turbulent Shear Flows 

William W. Liou 

1. Motivation and Objective 

This report documents the current progress in the research and development of 
modeling techniques for turbulent shear flows. These include a two-scale model for 
compressible turbulent flows and a new energy transfer model. The former rep- 
resents the status of our efforts to identify compressibility effects in turbulence. 
The energy transfer model refines a weakly nonlinear wave model developed ear- 
lier, which models directly the turbulent large structures. The objective of these 
activities is to develop second-order closures for compressible turbulent flows. 

2. Work Accomplished 

2.1 A Two-Scale Model for Compressible Turbulent Flows 

Numerical simulations of 2D and 3D compressible turbulence have shown that 
the existence of shocklet structures and the energy transfer mechanism between 
the kinetic energy and the thermal-energy are the two important compressibility 
effects 1,2,3,4 . These compressibility effects are incorporated into a new two-scale 
model. The model is based on the proposition that the effect of compressibility in 
turbulence is mainly on the energetic large eddies in turbulent shear flows. The small 
eddies are affected only indirectly through the increased spectral energy transfer. 
The development of the model and some results of its application to compressible 
free shear layers are briefly described here. A more detailed analysis is included in 
a NASA TM 5 . 

Firstly, it is assumed that the shocklet structures that may occur intermittently 
in compressible turbulent flows are formed mainly by the collision of the energetic 
turbulent eddies of large scale. The small eddies, which contain much less energy, 
are less efficient in the formation of shocklet structures when they collide with 
other eddies. Thus, the eddy shocklets scale with the energy containing eddies 
and have more direct influence on the evolution of the large eddies than on the 
smaller ones. The large vortical structures are intensified as they pass through the 
shocklet. This process, in other words, enhances the vortex stretching mechanism 
and increases the spectral energy transfer. In addition to the usual route of the 
vortex stretching mechanism that has already been enhanced, the small eddies may 
be generated directly from the passage of the large vortical structures through shock 
waves These processes of enhanced energy transfer may then cause the spectrum 
to depart from equilibrium. Another mechanism that may also contribute to the 
non-equilibrium spectrum or the creation of vorticity is strongly related to the 
pressure fluctuation. It has been shown by Kida and Orszag 3 and Lee et al. 2 , 
among others, that substantial vorticity is created by the baroclinic terms. The 
creation of vorticity, however, occurs mainly at the shock wave. 
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Based on the picture described above, the effect of compressibility in turbulence 
is mainly on the energy containing large eddies or the low wavenumber fluctuations. 
The large eddies respond more readily to changes in the compressible mean flow 
resulting from either high speed or combustion. The straining of the large eddies due 
to compressibility effects increases the spectral energy transfer to the small scales 
through the mechanism of vortex stretching and direct generation. The smaller 
scales, on the other hand, are only indirectly affected by compressibility. The energy 
contained in the small scales in the high wavenumber part of the energy spectrum 
is increased only as more energy is pumped in from the large eddies associated with 
the low wavenumber part of the spectrum. To model the Favre-averaged mean 
compressible turbulent quantities associated with these two distinct regimes in the 
energy spectrum we solve the modeled transport equations for the kinetic energy 
of the large eddy (k p ) and the small eddy (k t ) and the rate of energy transfer from 
the large eddy to the small eddy (e p ) and the rate of energy dissipation ( e t ). The 
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P.D. and E.S. denote the effects of pressure-dilatation and eddy shocklets, respec- 
tively. The definition of the model constants can be found in the NASA TM. The 
present two-scale model for compressible turbulence is built upon a parallel model 
for incompressible flows, Duncan et al. 6 . Models for the terms responsible for the 
compressibility effects are needed to close the equations. In this analysis, we have 
adopted Sarkar’s 7 model for the pressure-dilatation terms. To model the effects of 
the increased spectral energy transfer due to compressibility, a simple model has 
been constructed through dimensional reasoning. Its coefficient has a Mf depen- 
dence, similar to the dilatation dissipation model proposed by Zeman 8 and Sarkar 
et al. 9 . The compressibility corrections that they proposed have been implemented 
successfully into k - e models, Viegas and Rubesin 10 , into k-w models, Wilcox 11 
and into second-order closure models, Speziale and Sarkar 12 . 

Fig. 1 shows the variation of the vorticity thickness growth rate, d^/dx, as a 
function of convective Mach number. The vorticity thickness, <L, is defined by 

_ Uf-U, (5) 

w ~ (dU/dy)max K 1 
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The convective Mach number is defined as the ratio of the average convective 
velocity of the dominant large scale structures, relative to the free stream, to the free 
stream speeds of sound, Papamoschou and Roshko 13 . The convective Mach number 
has been shown to be an appropriate parameter to correlate experimental data and 
to identify the effects of compressibility. The vorticity thickness growth rates for 
compressible free shear flows, ( d8 u> /dx)(M c ,U s /Uf,p s /pf ), have been normalized 
by the corresponding values for incompressible flows, (d8u/dx) i {0,U s /Uf,p s /pf), 
and are presented in Fig. 1. The value of (d8 w /dx)i is obtained by using a relation, 
Papamoschou and Roshko 13 , 


d^ C 1 - ^■)( 1 + (^ 7 ) 1/2 ) 
d * 1 + ^(£) 1/2 


( 6 ) 


The constant of proportionality is obtained by the present model calculations per- 
formed in the limit of M c —* 0. Measured data are denoted by open symbols in Fig. 
1. Without the compressibility corrections, the current two-scale model and the 
two-scale model developed by Kim and Chen 14 (KC) predict a large reduction of 
the growth rate only at very high convective Mach numbers. With the inclusion of 
the effects of eddy shocklets and the pressure work, the current compressible two- 
scale model predicts a smooth reduction of the vorticity thickness growth rate as 
the convective Mach number increases. The calculated growth rate curve levels off 
at high convective Mach numbers. It should be noted that in the present analysis 
the convective Mach number of the shear layer is increased by increasing the Mach 
number of the high speed stream. According to the definition of the convective 
Mach number, there exists a maximum convective Mach number for a plane mixing 
layer of the same fluid with matched total temperature. That is, 
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where r = U s /Uf and 7 denotes the ratio of the specific heats of the working fluid. 
For a value of 12=0.1, the limiting convective Mach number for a plane shear layer 
of air is about 2.0. 


Since it is the Reynolds shear stress that appears in the mean momentum equa- 
tions and influences directly the development of the mean flow, it is interesting to 
see how its peak value varies as a function of M c . Note that in the current analysis, 
the Reynolds shear stress is related to the mean flow by a turbulent eddy viscosity, 
p t . That is, 
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In Fig. 2, the peak Reynolds shear stresses predicted by the present compressible 
two-scale model are compared with measured data, Elliott and Samimy 1 A The 
predictions of the present model, without compressibility corrections and of the 
KC model are also shown for comparison. The results show that, without the 
inclusion of some forms of compressibility corrections, none of the two-scale models 
tested, including the current model and the KC model, performs satisfactorily in the 
calculations of compressible free shear layers. The present compressible two-scale 
model under-predicts the absolute value of the peak Reynolds shear stress. However, 
the trend observed in the experiment that the level of the peak Reynolds shear stress 
decreases with increasing convective Mach number is picked up consistently by the 
current compressible two-scale model. Note that the Reynolds shear stress has 
been normalized by the square of the velocity difference of the two free streams. 
The model also shows that the value of the peak Reynolds shear stress appears to be 
independent of the velocity ratio of the free streams. In fact, the predicted variation 
of the peak value of the Reynolds shear stress as a function of the convective Mach 
number is similar to the predicted variation of the normalized vorticity thickness 
growth rate as a function of the convective Mach number. This characteristic of 
the present compressible two-scale model is consistent with the observation made 
by Elliott and Samimy 15 . They argue that, based on an integral analysis, the 
decreasing trend of the level of the Reynolds shear stress, as the convective Mach 
number is increased, is due mainly to the decrease of momentum thickness growth 
rate. However, the two speed ratios considered here, 0.1 and 0.2, are nearly equal 
to each other. Cases with a wider range of operating conditions, such as the speed 
ratios and the working fluids, need to be examined before any conclusive statement 
can be made. 

To further validate the present compressible two-scale model, it is applied to the 
compressible free shear layer corresponding to the Case 1 in Samimy and Elliott 16 . 
In this case, a fully expanded plane shear layer of air with M c = 0.51 and r = 0.36 is 
examined. The calculated mean profile shown in Fig. 3 agrees reasonably well with 
the measurement. As described previously, there are many possible causes for the 
small difference in the outer region of the mixing zone. Fig. 4 shows the comparison 
of the computed and the measured Reynolds shear stress. The present two-scale 
model under-predicts the peak Reynolds shear stress by about 12%. The profile of 
the Reynolds shear shear stress, however, agree very well with the measurement. 

2.2 A New Energy Transfer Model for Turbulent Free Shear Flows 

The model is built upon the weakly nonlinear wave models developed by Liou 
and Morris 17 . The development of the energy transfer model and some results of 
its application to an incompressible free shear layer are briefly described here. A 
more detailed analysis is included in a NASA TM 18 . 

The random flow properties are split into three components, 

fi = Ft + fi + (10) 

The fluctuation with respect to the long time-average component, Fi, is separated 
into a component representing the large-scale motion, fi, and one representing 
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the residual fluctuations, f The long time-average of the instantaneous value 
is denoted by an overbar, 



For the large-scale fluctuation, a separable form of solution is assumed, 

{u, v, p} = A(x) [u(y), v(y),p(y)] exp [i(ax - ut)\. (12) 


The bold face quantities denote a complex solution whose real part describes the 
physical properties of the large-scale structures, a (= a r + icti ) denotes a complex 
wavenumber and co the frequency. The governing equations for the local distribu- 
tions of the large structures can be reduced to the Rayleigh equation in terms of 

v, 

{(all - - a 2 ) - Q^} * = 0 (13) 

The amplitude, A(x), appears as a parameter in the local calculation for the u,v,p 
and is determined separately from the large scale turbulent kinetic energy equation, 
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where k = \uffui. k denotes the turbulent kinetic energy of the large-scale struc- 
ture. Note that in this analysis k denotes the turbulent kinetic energy of large 
scale structure of a single mode. The k p defined in the first part of this report 
represents the sum of the turbulent kinetic energy of all the modes in the entire 
large-scale spectrum. <> represents a short time- average with an average inter- 
val much smaller than T\ but much larger than the characteristic time scale of the 
background small-scale fluctuation, Strange and Crighton 19 . The interaction terms, 
the third expression on the right hand side of equation (14), describe the transfer 
of large-scale energy, presumably, to the small scales where energy is eventually 
dissipated by viscosity. The detailed analysis of the weakly nonlinear wave models 
and the numerical solution procedure used here can be found in Liou and Morris 17 . 


The spectral energy transfer results from the interactions between turbulent fluc- 
tuations of different scales. For the weakly nonlinear wave turbulence models, the 
energy transfer is of crucial importance in the determination of the wave ampli- 
tude and needs to be considered carefully. Very little information, experimental or 
theoretical, is available regarding the stresses, — < u[u'j >. 

The weakly nonlinear analysis seeks normal mode solution of the large-scale tur- 
bulent fluctuation. Locally, the fluctuations are described by the linearized Euler 
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equations. On the other hand, the spatial extent of each of the modes of the large- 
scale structures could be regarded as being determined by the wavenumber, ot r . 
Therefore, the proposition here is to estimate the characteristic size of the large 
scales as the wavelength associated with the large structures, which are predicted 
by the weakly nonlinear analysis. That is, 

l = l w = —. (15) 

Oi r 

where l w denotes the wavelength. Through dimensional reasoning, the enregy trans- 
fer can be modeled by 

jfcl 

C 2 f (16) 

This is the proposed model for the energy transfer from the large scale to the small 
scale. This estimate is in accord with the classic assumption of turbulence theory 
that dissipation “.. proceeds at a rate dictated by the inviscid inertia behavior of 
the large eddies”, Tennekes and Lumley 20 . Computationally, since the wavenumber 
is already a part of the solution of the equations for the large-scale fluctuation, this 
model involves no extra efforts in estimating the characteristic size of the energy 
containing large scales. This rather simple model provides a closure to the equations 
for the large-scale structure, thereby allowing render the weakly nonlinear wave 
description of the large-scale structure to be self-contained. This self-contained 
nature of the weakly nonlinear wave turbulence models may be important in the 
future applications to other turbulent free shear flows. 

The model is tested against an incompressible plane mixing layer. Since the most 
unstable mode interacts most strongly with the mean flow 17 , the most amplifying 
local instability is used in the modeling of the average, overall interactions between 
the mean and the large scale motions. Therefore, in the present formulation, the 
characteristic length scale l w is determined only by the locally most unstable modes. 

Fig. 5 shows the predicted evolution of the streamwise mean velocity profiles 
with axial distance, r] is a similarity coordinate, 


V ~ 2 / 1/2 

X — Xq 


(17) 


where jq / 2 denotes the location where the local mean velocity is one half of the free 
stream velocity. The predicted self-similar profiles agree well with that compiled 
by Patel 21 except at the low speed edge of the layer. Similar differences were also 
observed by Liou and Morris 17 . They attributed this difference to the single mode 
representation of the entire large scale spectrum and the uncertainties in the mea- 
surements in this region resulting from the local large changes in the instantaneous 
flow direction. 

The streamwise evolution of the amplitude of the large-scale structures is shown 
in Fig. 6. After a region of establishment, the amplitude reaches a saturated value. 
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In this region, the rate of the production of the large-scale turbulent kinetic energy 
from the mean flow is balanced by the rate of energy transfer from the large scales 
to the small scales. Note that, for the present energy transfer model, the amplitude 
equation becomes, 

^ = G 3 (x) A 2 - G 4 (x) A 3 (18) 

ax 

G 3 and G4 denote the normalized positive definite integrals of the production terms 
and interaction terms across the layer, respectively. The critical points of the nonlin- 
ear equation (18), where dA 2 /dx = 0, are A\ = 0 and G 4 (x 2 )A 2 = G 3 (x 2 )- Simple 
analyses by applying the Liapunov function method 22 show that A\ is an unstable 
critical point. Any small disturbances to A u say A[ would grow exponentially. In 
fact, 

(A 2 )' « e G3(ll) 1 (19) 

A 2 , on the other hand, is asymptotically stable. A disturbance about the A 2 , say 
A' 2 , would decay exponentially, 
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The saturated value of the amplitude, A 2 , is an asymptotically equilibrium value. 
It indicates an asymptotically equilibrium state of the large-scale structures. The 
simple instability analyses also show that any deviation away from this equilibrium 
state would be damped out exponentially. Consequently, the saturation of the wave 
amplitude may provide an indication of the the self-similarity of the flow in terms 
of the development of the large-scale structures. 


3. Future Plans 

3.1 A Two-Scale Model for Compressible Turbulent Flows 

(1) Extend the two-scale model to wall-bounded flows. 

(2) Continue the development of second-order closure models that account explic- 
itly for the compressibility effects identified during the development of the two-scale 
eddy- viscosity model. 

3.2 A New Energy Transfer Model for Turbulent Free Shear Flows 

(1) Apply the weakly nonlinear wave model to compressible mixing layers to 
investigate the effects of compressibility on the characteristics of the coherent large- 
scale structures. 
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Figure 1. Variation of relative growth rate with convective Mach number, r = 0.1. — , 

Present: without compressibility corrections; , Present: with compressibility 

corrections; , KC; ° , Papamoschou and Roshko; A , Goebel et al ; <> , 

Samimy and Elliott; , Sullins et al , , Messersmith et al ; X , Chinzei et al ; , 

Ikawa and Kubota et al. 
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Figure 2. Variation of the peak Reynolds shear stress with convective Mach number. — £ — 
Present: r = 0.1, without compressibility corrections; — S — , Present: r = 0. 

with compressibility corrections; *• — , Present: r = 0.2, without compressibilit 

corrections; 1 , Present: r = 0.2, with compressibility corrections; , K( 

r = 0.1; • , Elliott and Samimy. 
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Normalized mean velocity, U*. 
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Figure 3. Comparison of mean velocity profiles for the shear layer in Case 1 of Elliott and 
Samimy (1991). , Present; • , experiment. 
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Figure 4. Comparison of Reynolds shear stress for the shear layer in Case 1 of Elliott and 
Samimy (1991). , Present; • , experiment. 
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Figure 5. Mean velocity profiles. , x = 3.82; , 4.43; , 5.99; 


0.20 


• , Patel 9 . 
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Figure 6. Variation of the wave amplitude with streamwise distance. 
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Modeling of Near Wall Turbulence 
and Modeling of Bypass Transition 

Z. Yang 

1. Motivation and Objective 

1) Modeling of the near wall turbulence. We aim to develop a second order clo- 
sure for near wall turbulence. As a first step of this project, we try to develop a 
k - e model for near wall turbulence. We require the resulting model to be able to 
handle both near wall turbulence and turbulent flows away from the wall, compu- 
tationally robust, and applicable for complex flow situations, flow with separation, 
for example. 

2) Modeling of the bypass transition. We aim to develop a bypass transition 
model which contains the effect of intermittency. Thus, the model can be used 
for both the transitional boundary layers and the turbulent boundary layers. We 
require the resulting model to give a good prediction of momentum and heat transfer 
within the transitional boundary and a good prediction of the effect of freestream 
turbulence on transitional boundary layers. 

2. Work Accomplished 

In the past year, progress has been made in both topics mentioned above (i.e., 
modeling of the near wall turbulence and modeling of the bypass transition). In the 
paragraphs below, these two topics will be reported separately. 

2.1 Modeling of Near Wall Turbulence 

Because of the wide range of scales involved in a turbulent flow, DNS (direct 
numerical simulation) is limited to flows of moderate Reynolds number and simple 
geometry. Turbulence modeling is the only viable approach for the calculation of 
turbulent flows of engineering interest. In turbulence modeling, the k-e model is the 
most widely used model in engineering calculations. The Standard k — e Model 12 
was devised for high Reynolds number turbulent flows and is traditionally used 
in conjunction with a wall function when it is applied to wall bounded turbulent 
flows. However universal wall functions do not exist in complex flows and it is thus 
necessary to develop a form of k — e model equations which can be integrated down 
to the wall. 

Jones and Launder 3 were the first to propose a low Reynolds number k-e model 
for near wall turbulence, which was then followed by a number of similar k — e 
models. A critical evaluation of the pre-1985 models was made by Patel et al. 4 . More 
recently proposed models can be found in Shih 5 and Lang and Shih 6 . Three major 
deficiencies can be pointed out about existing k-e models. (Some of the models may 
have only one or two of the three deficiencies.) First, a near wall pseudo-dissipation 
rate was introduced to remove the singularity in the dissipation equation at the 
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wall. The definition of the near wall pseudo-dissipation rate was quite arbitrary. 
Second, the model constants were different from those of the Standard k — e Model, 
making the near wall models less capable of handling flows containing both high 
Reynolds number turbulence and near wall turbulence, which often occurs for real 
flow situations. Patel et al. 4 put as the first criterion the ability of the near wall 
models to predict turbulent free shear flows. Third, the variable y + is used in the 
damping function / M of the eddy viscosity formulae. Since the definition of y + 
involves u T , the friction velocity, any model containing y + cannot be used in flows 
with separation. 

Effort is made to propose a new k-e model for near turbulence which is free of the 
three deficiencies mentioned above. In this model, k 1//2 is chosen as the turbulent 
velocity scale. The time scale is bounded from below by the Kolmogorov time scale. 
The dissipation equation is reformulated using this time scale and no singularity 
exists at the wall. Thus, it is no longer necessary to introduce the near wall pseudo- 
dissipation rate. The model constants used are the same as in the Standard k-e 
Model. Thus, the proposed model will be also suitable for flows away from the wall. 
An earlier version of the model, which contains y + in the damping function, was 
proposed and reported in Yang and Shih 7 . The model is now improved by using 
R y = k'^yjv instead of y + in the damping function. Hence, the present model can 
be used for flows with separation. 

In the present model, the eddy viscosity is given by 


vt — c^f^kT ( 1 ) 

where the time scale T is written as 

( 2 ) 

The first part is the time scale conventionally used for high Reynolds number tur- 
bulent flows and the second part is the Kolmogorov time scale. Away from the wall, 
the first part is much large than the second part while near the wall the second part 
dominates, giving the Kolmogorov time scale as the turbulent time scale at the wall. 
The time scale given is bounded from below by the Kolmogorov time scale and is 
always positive. 

The damping function is given by 


/ M = [1 — exp{—aiR y — a$Ry — a 5 -Ry )] 1 ^ 2 (3) 

where a x = 1.5 x 10“ 4 , a 3 = 5.0 x 10~ 7 , a 5 = 1.0 x 10" 10 . The damping function 
is chosen such that the shear stress has the correct near wall asymptotic behavior. 
Away from the wall, approaches one as required. 

The modeled transport equations for k and e are 

k + Ujkj = [( 1 / + — )k,j],j- < UiUj > Uij - e, (4) 
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e + Ujtj = [(u + — )e,j]j + (— C\ e < u-iUj > Uij - Cie^fT + vvTUi,jkUi,jk- (5) 

^ i 


Because T is always positive, the dissipation rate equation does not have a singu- 
larity at the wall. 

Turbulent channel flows and flat plate boundary layer flows at different Reynolds 
numbers were calculated using the proposed model. The mean velocity, turbu- 
lent kinetic energy, turbulent shear stress and turbulent dissipation rate for tur- 
bulent channel flow at Re r = 395 and turbulent flat plate boundary layer flow at 
Reo = 1410 are shown in Fig 1 and Fig 2, respectively. DNS data for these cases 
are shown for comparison. Also shown are the predictions using the Jones-Launder 
model and the k - e model proposed by Chien 8 . These two models are chosen 
because the Jones-Launder model is the first k - e model for near wall turbulence 
while Chien’s model is known to perform quite well for turbulent boundary layer 
flows. Overall, the proposed model is found to give a better prediction. Calcula- 
tions were also made for turbulent flat plate boundary layers at larger Reynolds 
numbers, turbulent boundary layers with pressure gradient (favorable pressure gra- 
dient, adverse pressure gradient, and increasingly adverse pressure.) The results of 
these computation and the comparisons with the available experimental data can 
be found in Yang and Shih 9 . 


2.2 Modeling of Bypass Transition 

In a quiescent environment, transition is preceeded by the amplification of Tollmien- 
Schlichting waves. These waves eventually break down, giving rise to turbulent 
spots, which can be viewed as the onset of transition. In an environment with high 
freestream turbulence, say the flow passing over a turbine blade, turbulent spots 
are formed due to the transport of turbulence from the freestream to the boundary 
layer rather than the T-S wave amplification. This type of transition is called bypass 
transition. Accurate prediction of bypass transitional boundary layers is very im- 
portant for internal fluid mechanics because a significant proportion of the turbine 
blade is in the a transitional boundary layer region. Furthermore, the performance 
and the life of a turbine are directly related to the peak vales of the momentum and 
heat transfer both of which occur in the transitional boundary layer. 

Priddin 10 was the first to notice that the low Reynolds number two equation 
models have the potential to predict transitional flows under the influence of the 
freestream turbulence. This is probably due to the fact that the generation of 
turbulent spots in a boundary layer is a random process and the flow is almost 
fully developed turbulent within a turbulent spot. A detailed calculation procedure 
was given by Rodi and Scheuerer 11 , in which the Lam & Bremhorst low Reynolds 
number k - e model was used. More recently, a comparative study of the perfor- 
mance of existing low Reynolds number k — e models in predicting laminar-turbulent 
transition was made by Fujisawa 12 . 

While the low Reynolds number k - e models could mimic transition, the quan- 
titative predictions do not compare very well with the experimental data. This is 
due to the fact that all these low Reynolds number k - e models were originally 
proposed for fully developed turbulent flows and did not take into consideration the 
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distinct feature of a transitional boundary layer - intermittency. The intermittency 
of a transitional boundary layer is measured by the intermittency factor which can 
be viewed as the percentage time a transitional boundary layer is in the turbulent 
state due to the passing of a turbulent spot. 

We propose a model for the calculation of transitional boundary layers, which 
takes the effect of intermittency into consideration. The model is based on the k — c 
model for near wall turbulence we have stated above. The effect of intermittency 
is introduced through the following argument: since the percentage of time that a 
transitional boundary layer is turbulent is measured by the intermittency factor, 
a model for transitional boundary layers could be constructed from a model for 
turbulent boundary layers by multiplying all the terms due to turbulence mechanism 
by a weighting factor 7, which is linearly related to the intermittency factor. Thus, 
the governing equations for the flat plate transitional boundary layers are, after 
using the boundary layer approximation and the eddy viscosity assumption, 
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dx dy 
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where x, y are the coordinates along and normal to the plate and U, V are the mean 
velocities in the x, y directions, respectively. 

In order to close the above equations, an expression for the weighting factor 7 
is needed. The weighting factor is assumed to be related to both the freestream 
turbulent level and the intermittency factor of the boundary layer. The intermit- 
tency factor is assumed to be determined by the local state of the boundary layer. 
We use H, the shape factor, to characterize the local state of the boundary layer 
since both the intermittency factor and the shape factor change monotonically from 
the laminar boundary layer to the turbulent boundary layer. Experimental results 
by Abu-Ghannam and Shaw 13 are used as a guide to construct this function. The 
weighting factor is also assumed to change in the y direction in such a way that 
outside the boundary layer, the weighting factor is one since the freestream is gov- 
erned by decaying turbulence. The final expression for the weighting factor is given 
in Yang and Shih 14 . 

The above system of parabolic equations need to be supplymented by boundary 
conditions at the wall and at the freestream and by initial conditions at the starting 
point of the calculation. At the wall 


U = V = k = 0, 
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The boundary condition for e was obtained by applying the equation for turbulent 
kinetic energy down to the wall. 

At the edge of the boundary layer, the flow variables are given by the freestream 
values, i.e. 

U = U e ,k = k e ,e = e e . (11) 


The k e ,e e in the above are found from the transport equations for k and e, with 
the condition that the gradient of the flow variables in the y direction vanishes as 
the free stream is approached. Thus, 


dk P 

U e - r= 

ax 
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rfe. _ Cute (13) 

e dx T 

k e o and e e0 (the values of k e and c e at the leading edge, for example) are needed. 
k e o is obtained from the experiment, and e e0 is determined in such a way that the 
resulting k e (x ) profile agrees with the experiment. 

One of the issues in the calculation of transitional boundary layers through the 
low Reynolds number k - e models is the prescription of the initial profiles for 
the turbulence kinetic energy and its dissipation rate, the later of which could not 
be found from the experiment directly. An expression for the initial profiles were 
given in Rodi and Scheuerer 11 . However, computations by Yang and Shih 15 which 
tested the effect of the initial conditions on the transition prediction found, in 
agreement with the findings of Patankar and Schmit 16 , that the predicted onset of 
the transition is sensitive to the initial profiles. This sensitivity of the results to the 
initial conditions suggests that the only place where the initial conditions could be 
specified unambiguously is at the leading edge. At the leading edge, tne turbulent 
kinetic energy and its dissipation rate take constant profiles, the values of which aie 
determined by the law for the decaying turbulence. 

With the initial conditions given at the leading edge and the boundary conditions 
given above, the solutions are marched downstream. Flat plate boundary layers with 
free stream turbulence levels of 3% (Case T3A) and 6% (Case T3B) respectively 
were calculated using the present model. These are the benchmark cases in an 
ongoing project coordinated by Savill 1 ' , testing the capability of turbulence models 
in predicting transitional flows. Fig. 3 shows the variation of skin friction coefficient 
C f against Re x . Results from the experiment is shown for comparison. In addition, 
the prediction of the Launder-Sharma model is also shown in the figure because it 
was reported that among the lower Reynolds number k - e models, the Launder- 
Sharma model performs best for transitional boundary layers. It is clear that the 
present model gives a better prediction. Other features in the transitional boundary 
layers and the calculations of the transitional boundary layers with other levels of 
freestream turbulence can be found in Ref. 14. 
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3. Future Plans 

1) Modeling of the near wall turbulence. 

The proposed k - e model is only tested for simple parabolic flows so far. Because 
of the form of the model equations, the proposed model can be used in complex 
flow situations, flow with separation for example. The performance of the proposed 
model in those situations will be tested. 

We will work on the second order closure for near wall turbulence. In particular, 
we will be looking at the effect of the mean flow inhomogeneity on the pressure 
strain correlation. We are hoping to represent this effect rationally, so that the 
ad hoc damping functions currently being used in all the near wall second order 
closures can be avoided. 

2) Modeling of bypass transition. 

We will apply the proposed model to transitional boundary layers with pressure 
gradient and curvature. We will also extend the model to thermal boundary layers. 
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Fig la: Mean velocity profile for channel flow at Re r = 395. 
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Fig lb: Turbulent energy profile for channel flow at Re r = 395. 
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Fig lc: Shear stress profile for channel flow at Re T — 395. 



0 0 25.0 50.0 75.0 100.0 

y + 

Fig Id: Dissipation rate profile for channel flow at Re r = 395. 
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Fig 2a: Mean velocity profile for channel flow at Re$ = 1410. 
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Fig 2b: Turbulent energy profile for channel flow at Re$ = 1410. 
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Fig 2c: Shear stress profile for channel flow at Re$ = 1410. 



Fig 2d: Dissipation rate profile for channel flow at Reg = 1410. 
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Fig 3a: Skin friction variation for T3A. 



Fig 3b: Skin friction variation for T3B. 
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Two Equation Modelling and the 
Pseudo Compressibility Technique 

C. J. Steffen, Jr. 

1. Motivation and Objective 

The primary objective of the Center for Modelling of Turbulence and Transition 
(CMOTT) is to further the understanding of turbulence theory for engineering 
applications. One important foundation for this research is the establishment of 
a data base encompassing the multitude of existing models as well as newly proposed 
ideas. The research effort described in the next few pages is a precursor to an extended 
survey of two equation turbulence models in the presence of a separated shear layer. 

Recently, several authors have examined the performance of two equation models 
in the context of the backward facing step flow. Conflicting results, however, demand 
that further attention is necessary to properly understand the behavior and limita- 
tions of this popular technique, especially the low Reynolds number formulations. 
The objective of this research is to validate an incompressible Navier Stokes code 
for use as a numerical test-bed. In turn, this code will be used for analyzing the 
performance of several two equation models. 

2. Work Accomplished 

To date, the validation of the incompressible code DTNS is complete. The details 
of this validation study are documented in reference[l]. The code is based upon 
the pseudo-compressibility technique and incorporates the approximate factorization 
scheme for time integration. Two laminar benchmark flows are used to measuie the 
performance and implementation of the numerical methods. The classic solution of 
the Blasius boundary layer is used for validating the flat plate flow, while experimental 
data is incorporated in the validation of backward facing step flow. 

An initial result for the standard high Reynolds number k-e equations has also been 
calculated to demonstrate an initial performance level of the solution technique for 
the turbulence equations. 

2.1 Numerical Method 

The absence of a variable density in the continuity equation governing incompress- 
ible flow complicates the numerical integration procedure. One solution technique 
that has been well received is that of pseudo-compressibility. This idea was first 
put forward by Chorin[2] and enables the equations to be solved using the prim- 
itive variables. Recently, Chang and I\wak[3], Rizzi and Eiiksson[4], Kwak and 
Chakravarthy[5], Michelassi and Shih[6], and Turkel[7] have found this method suit- 
able for resolving incompressible flow. This particular implementation has been vali- 
dated by Gorski[S-10] for several different benchmark flows. 
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Examine the system of equations involved in the pseudo compressibility method 
and notice that they differ from the incompressible Navier Stokes equations by the 
addition of a time dependant pressure term in the continuity equation. 
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Here, x and y are the independent variables and Re refers to the Reynolds number. 
The constant (3 is known as the pseudo compressibility parameter. This system 
is hyperbolic in nature whi le the incompressible flow equations are elliptic. The 
pseudo sound speed, c = yj u 2 + (3, is governed by the value of the parameter /?, 
whereas the physical sound speed is infinite. Chang and Kwak [3] have shown that 
for (3 > 0 the finite speed pseudo waves vanish as time progresses and yield the proper 
incompressible solution at the steady state limit. It is through this parameter /3 that 
the convective and acoustic waves are decoupled, and thus convergence is governed. In 
choosing an optimum value for this parameter, the goal is to avoid giving the viscous 
effects time to react to the passing of the nonphysical transient pressure waves. Thus 
a lower bound on the acoustic speeds translates into a lower bound on f3. However, 
an upper bound on (3 is strictly scheme dependent. 

The approximate factorization is rather straight forward. Equations (2) can be 
rewritten as 
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or in generalized coordinates (£(x,y),T}(x,y)) as seen here in equation (3). 
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The third order accurate Chakravarthy and Osher TVD scheme is used to discretize 
the convective terms while central differencing is used for the viscous fluxes. Thus 
the approximate factorization of this method is expressed in terms of flux Jacobians 
A and B as 
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where V is the cell volume, and At = is the timestep. The matrices A„ and B„ 

are the viscous flux Jacobians. The interested reader is encouraged to see reference^] 
for more details. 

The k-e equations are solved in a similar manner as that seen in equation 2. 
Consider the standard high Reynolds number form of the transport equations of the 
turbulent kinetic energy and dissipation of turbulent genetic energy as in equation 
(5): 
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and the production of turbulent genetic energy defined as shown below. 
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Again, we can translate these equations in a manner similar to that of equation (3) 
and create the following system for generalized coordinates: 
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This can be solved via approximate factorization as shown here in equation 9. 
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Again, the scheme is written in terms of the turbulence flux Jacobians and the 
approximate source Jacobian, H*\ The coefficients a ^ and a ^ are provided to enhance 
the diagonal dominance of the implicit method. While the convective flux terms are 
linear, the linearization of the source terms is not straight forward. Furthermore, the 
low Reynolds number formulation of equation 5 includes other nonlinear correction 
terms in the source vector S l . The interested reader is encouraged to see the recent 
paper by Michelassi and Shih[6]. 


2.2 Discussion 

Overall the performance for laminar flow documented in reference[l] is very en- 
couraging. The flat plate results shown agree with the theoretical Blasius solution 
for several different grid configurations. The laminar back facing step flow was also 
well resolved both in terms of the primary recirculation zone reattachment length, 
and corresponding velocity profiles as seen in figures 1 and 2. However, performance 
of the two dimensional flow code drops off dramatically when the Reynolds number 
corresponding to the onset of 3D flow is exceeded. Armaly et al [11] observed this 
phenomenon occur at a Reynolds number of 400, based upon the hydraulic diameter 
of the inlet channel. Again, the numerical results bear this finding out. 


Figure 1 Reattachment length predictions — comparison with experiment [11], xi, X 2 , and X3 are the 
primary reattachment, secondary detachment, and secondary reattachment lengths, respectively. 
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Figure 2 ReD=389 velocity profiles — comparison with experiment [11]. 
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The current thrust of effort is in establishing a comparison between the standard 
high Reynolds number form k-e model and a version which incorporates the low 
Reynolds number corrections, such as that of Jones and Launder[12]. This model 
is particularly convenient in that the damping functions are independent of the 
parameter, y + . Thus, multiple walls bounding the computational domain do not need 
special attention. An initial calculation using the high Reynolds number approach 
agrees with the published results of Speziale and Thangam[13]. However, the Jones 
and Launder model is proving to be quite challenging. 

3. Future Plans 

Eventually, several two equation models will be examined for their behavior in the 
back facing step flow. The following is a partial list of the two equation models of 
primary interest: 

1. Jones and Launder 

2. Chien 

3. Yang and Shih 

4. Shih and Lumley 

The current research effort will result in an increased understanding of the need 
for proper modelling of the near wall region. However, the increase in accuracy 
maybe offset by an unwarranted increase in effort. Further numerical experiments 
are necessary before any conclusions can be drawn. 
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Renormalization Group Methods for the 
Reynolds Stress Transport Equations 

R. Rubinstein 


1. Motivation and Objective 

The Yakhot-Orszag renormalization group is used to analyze the pressure gradient- 
velocity correlation and return to isotropy terms in the Reynolds stress transport 
equations. The perturbation series for the relevant correlations, evaluated to lowest 
order in the e-expansion of the Yakhot-Orszag theory, are infinite series in tensor 
product powers of the mean velocity gradient and its transpose. Formal lowest 
order Pade approximations to the sums of these series produce a rapid pressure 
strain model of the form proposed by Launder, Reece, and Rodi, and a return to 
isotropy model of the form proposed by Rotta. In both cases, the model constants 
are computed theoretically. The predicted Reynolds stress ratios in simple shear 
flows are evaluated and compared with experimental data. The possibility is dis- 
cussed of deriving higher order nonlinear models by approximating the sums more 
accurately. 

The Yakhot-Orszag renormalization group provides a systematic procedure for de- 
riving turbulence models. Typical applications have included theoretical derivation 
of the universal constants of isotropic turbulence theory, such as the Kolmogorov 
constant, and derivation of two equation models, again with theoretically com- 
puted constants and low Reynolds number forms of the equations. Recent work 
has applied this formalism to Reynolds stress modeling, previously in the form of a 
nonlinear eddy viscosity representation of the Reynolds stresses, which can be used 
to model the simplest normal stress effects. The present work attempts to apply 
the Yakhot-Orszag formalism to Reynolds stress transport modeling. 

2. Work Accomplished 

The modelling of the pressure gradient-velocity correlation and return to isotropy 

term in the Reynolds stress transport equation remains an area of active research . 1 ’ 2,3 

Models will be derived here using the Yakhot-Orszag renormalization group 2 * 4 par- 

tially along the lines of our previous work 5 . The result is a model for the rapid 
pressure-strain term of the form proposed by Launder, Reece and Rodi 6 (LRR) and 

a model for return to isotropy of the form proposed by Rotta 7 with theoretically 
computed constants in good agreement with accepted values. As is usual in inves- 
tigations of this sort, the priority of Yoshizawa in deriving a pressure strain model 
analytically 8 must be noted. 

The analysis requires some new ideas in renormalization group theory recently 
introduced by Yakhot et al 9 . As Yakhot et al 9 emphasize, the application of the 
renormalization group mode elimination formalism to shear flow creates a double 

perturbation series in powers of e, the parameter of the isotropic theory, and in 
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powers of a dimensionless strain rate, 77 = SK/e, where K denotes the turbulence 
kinetic energy, e denotes the dissipation rate, and S is a measure of the mean 
strain: in Ref. 9, S 2 = + §^ JL )^ aL - The present analysis also leads to double 

expansions of this type, with the powers S n replaced by tensors homogeneous of 
degree n in the mean velocity gradient VU and its transpose. It will be convenient to 
retain the terminology of Ref. 9 and call this expansion the 77-expansion; when the 
distinction is pertinent, the expansion of Ref. 9 will be called a scalar 77-expansion. 

The heuristic program of evaluating all scalar amplitudes to lowest order in e 
has proven successful in the past: apparently, the e-expansion is an asymptotic 
series with a sum given very nearly by its first term 10 . Unfortunately, there is no 
analogous basis for truncating the 77-expansion. There are fundamental reasons for 
this distinction between these expansions. The present 77-expansion is tensorial: 
successively higher order terms do not introduce merely numerical corrections, but 
increasingly complex asymmetries into the theory. Truncation therefore imposes 
a possibly inappropriate symmetry or other constraint on the model. Thus, in 
Ref. 5 the 77-expansion for the Reynolds stress r was truncated at second order as 
suggested by previous work of Yoshizawa 11 and Speziale 12 . Although this type of 
modelling permits unequal normal stresses in a simple shear flow, it is not maximally 
asymmetric: for example, in a flow with mean velocity components Uj(xi,x 2), a 
cubic model including a term r ~ VU 2 VU T + VUVt / T2 would permit nonzero t 23 
in the presence of vanishing dU^/dx^ and dU^/dx?., an effect which cannot be ruled 
out in advance. 

Although generalizations 13 of the Cayley-Hamilton Theorem limit the number 
of independent tensors S^ n ^, anisotropy and asymmetry cannot exist at all without 
some terms of higher order in 77; indeed, truncation at lowest order in 77 just produces 
a theory of isotropic turbulence. But the series truncated at any higher order can 
be unsatisfactory in flow regions in which some components (^)K/e are large. In 
such regions, the truncated series is dominated by its highest order terms. For the 
quadratic stress models of Refs. 5, 11, 12, this domination can produce negative 
normal stresses in the buffer layers of wall bounded flows. Increasing the order of 
truncation obviously exacerbates this problem. 

It follows that finite truncation of the 77-expansion is theoretically unsatisfactory. 
Yakhot et al 9 therefore propose that this expansion must be summed, even if only 
approximately, and have suggested a prototype summation in a different context. 
It should be noted that the same issues arise naturally in Yoshizawa’s formalism, 
which also generates infinite series in the mean velocity gradients (and in other 
quantities as well) for correlations of interest in turbulence modeling. Yoshizawa 
has concluded independently that summation of this series is essential and has also 
derived a Reynolds stress transport model by introducing such summations 14 . 

In this paper, the perturbation series which the Yakhot-Orszag renormalization 
group generates for the correlation 
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where (0) denotes deviatoric part and ( ij ) denotes index interchange in the preced- 
ing term. 

The next order will produce a term T3 containing cubic products of velocities u < . 
In view of the form of the LRR model, it is reasonable to ask whether a term with 
only one gradient, proportional in the high Reynolds number limit to tVU might 
occur at this order. Such terms do occur, but they cancel. Evaluation of T3 proves 
to require expansions of the projection operators to second order, leading instead 
to terms homogeneous of degree three in the mean velocity gradient and its 
transpose. In general, the term T n of order n has the form S ^ n \K/e) n . As noted 
in the Introduction, it will be imperative to include effects of all orders in SK/e 
in the model, but because the terms T n involve ever higher order derivatives of 
the transverse projection operators, they do not have an obvious law of formation. 
Therefore, an exact summation does not appear feasible. 

A simple approximate summation is obtained by introducing the perturbation 
series 5 for v^v~^ in the form 

and dropping the quadratic terms. The resulting model, 


n « = 5* 


dUi dUi 


dx, 


+ 


dxi 


4* Cr 


T 1 


UiUp 


( 0 ) 


+ cj -2 


'U'i'U'p 


(0)dUp 

dxj 


,0) dU p 

+ UjUpi dx. 


dUj 

dx v +UjUp 
-i(o) 


( 0 ) 


dUj 

dx 


1 ( 0 ) 


v J 


( 5 ) 


with 

Cri = ^ = -7619 Cr 2 = 1; = -0952 (6) 

agrees with the perturbation series (3) and (4) to terms of order S^ 3 ). However, 
unlike the explicit quadratic model which results from simply dropping the 0(S^ 3 ^) 
terms, this model includes effects of all order in SK/e. The consequences of this 
fact will be discussed later. This type of summation has also been applied by 
Yoshizawa 14 . Eqs. (5) and (6) can be compared with Launder, Reece and Rodi’s 
“model 1”, Eq. (5) with the empirically adjusted constants 


Cr 1 = .7636 Cr 2 = .1091 (7) 

In this model, the constants Cr 1 and C+ 2 were not chosen independently; instead, 
to insure some consistency conditions introduced by Rotta 7 , Launder, Reece, and 
Rodi set 6 


Cr 1 


C 2 + 8 
_____ 


Cr2 = 


8C2-2 
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where only the constant C 2 is arbitrary. By eliminating C 2 between these equations, 
there results 


8Cri -Cf 2 = 6 (8) 

which is also satisfied by the choice of constants in Eq. (6). The LRR model 
corresponds to the choice C 2 = .4; Eqs. (5) and (6) correspond instead to the 
choice C 2 = 8/21 ~ .36. 

The approximate summation used to derive Eq. (5) can be systematically gener- 
alized to generate an infinite number of models for IJij. For example, suppose that 
the perturbation series for r is introduced into the cubic terms in the perturbation 
series instead of in the quadratic terms as above. This substitution will produce a 
model which can be written symbolically in the form 

n ~ S (1) + S (2) + r(S (1) ' + S (2) ') 

where tS W' denotes a sum of matrix products in all possible orders of r and terms 
S (i) . The requirement that the original series agree to order S (4) with the approxi- 
mation when r is replaced by its perturbation series determines this approximation 
uniquely. 

2.2 The Return to Isotropy Model 

The analytical description of return to isotropy is no less controversial than the 
modeling of the fast pressure strain term 3 . In the usual approach to turbulence 
modeling, in which correlations generated by Reynolds averaging are closed phe- 
nomenologically, this process is considered to result partly from the pressure cor- 
relation through a “slow” term independent of the mean flow, and partly from the 
deviatoric part of the dissipative correlation ^of Fro™ tllis viewpoint, the 
analysis in Sect. I is incomplete because it discloses only a term proportional to the 
mean velocity gradient, but no slow term. The return to isotropy will be derived 
here by renormalization group methods following a suggestion of Yakhot 17 . 

From the renormalization group viewpoint, it is natural to investigate the return 
to isotropy, even independently of the stress transport equation, by writing the 
perturbation series for 

Ui lH +Uj lit = J Ui ^ ~ 4)(-iu)uj(q)dq + ( ij ) (9) 

This perturbation series differs from the perturbation series for the Reynolds stresses 
previously reported 5 only in the occurence of an additional factor —iu) in all fre- 
quency integrals. 

The analysis is straightforward. Only the deviatoric terms require attention be- 
cause the part of the correlation proporational to 6ij contributes to the transport 
equation for K which has been analyzed by Yakhot and Smith 15 . The lowest order 
deviator appears at first order in 77 ; to lowest order in e 
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where 


Tx 


(dUi dU£\ 

V dxj dxi / 


( 10 ) 


dt i _ 1 V 
dr 15 uA 2 


( 11 ) 


In View of the form of the Rotta model, it is reasonable to seek terms at the next 
order proportional to U{Uj. As in Sect. I, such terms do appear, but cancel exactly. 
This apparently ubiquitous cancellation was also obtained by Smith and Reynolds 18 
in an analysis of the e transport equation. Accordingly, the second order analysis 
in t] produces quadratic terms in the velocity gradients. Finite truncation of this 
series violates the requirement that return to isotropy be independent of the mean 
flow. Therefore, we must seek a reasonable approximate summation. The form of 
the lowest order term given in Eqs. (10) and (11) suggests 


n 'ii = 


tk v (Wi + d _!> 

v V dxj dx 


j\ 1 i 

— ~ UiU 

)i J V 


:( 0 ) 


Despite its triviality, this replacement does produce an approximate sum which 
agrees exactly with perturbation theory to lowest order. It therefore can be consid- 
ered a type of Pade approximation. 

At the infinite Reynolds number asymptotic limit 


n[j = -Cr^uTutW (12) 

where, in the Yakhot-Orszag theory, Cr = V/e ~ 1.6. Equation (12) is therefore 
simply the standard Rotta model with Rotta constant ~ 1.6 in agreement with an 
earlier proposal of Yakhot 17 . 

A preliminary discussion of higher order summation may be appropriate. By 
analyzing the spectral dynamics of the return to isotropy, Weinstock 3 concluded 
that the shear and normal stresses relax at different rates. Although this behavior 
is obviously not accommodated by the Rotta model, it is consistent with the present 
theory: the perturbation series for JJ' is obtained from the series for r by multiplying 
the term of order n by the factor C n £/K for some constant C n . The C n are all 
unequal; therefore, the Rotta model is not exact. Now comparison with the series 
for t shows 5 that relaxation of the shear stress is governed by the linear term 
whereas relaxation of the normal stresses is governed by the quadratic term S^ 2 ). 
Since C 2 7 ^ C\ , these stresses relax at different rates. The difference is suppressed 
in the Rotta model, which arose in the present formalism by replacing all of the C n 
by Ci. 

2.3 Algebraic Reynolds Stress Models 

The approximation, due to Rodi 20 , of the Reynolds stress transport equation 
by an algebraic model under the conditions of semi-homogeneous flow (negligible 
diffusion of r and t/K approximately constant) takes the form 
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P - e 

K 


7(0) _ _ 


UiUp ' = 



dUj 

dx p 


I UjUp 


dUj 

dxp 


( 0 ) 

+ n'ij + n tJ 


(13) 


where 77 and 77' depend on r and V77. Explicit solutions for r can be obtained, 
at least in principle, for any such approximation 25 . Briefly, one introduces a basis 
for polynomials in Vt/, and VC/ T . The basis contains 11 terms of homogeneity 
order n < 5. Writing r as a sum of these terms with unknown coefficients and 
substituting in Eq. (13) leads to the explicit expression 


t/K = Y H\ m) S\ n) (W, W T ) 


(14) 


where H^ n ' > is a scalar function of VU and VL rT such that 


77} m) ~ \VU\ m 

when | VC/| — > oo. The assumptions made on the approximate summations require 
m + n = 0; thus, r/K is bounded when SK/e — > oo. For example, the familiar eddy 
viscosity formula is replaced in Eq. (29) by a term 


(V(7, WU t ) (VF + V£/ t ) 

Pope observed 25 that the coefficients 77( _n ) in Eq. (14) would certainly be in- 
tractably complex; although they could be explicitly exhibited by symbolic com- 
putation, the result would only pertain to the particular implicit equation for the 
Reynolds stresses assumed initially in Eq. (13). Therefore, it is equally reason- 
able just to postulate simple forms for the functions 77( _n L This type of modeling 
could be particularly interesting when applied to the coefficients of the quadratically 
nonlinear models of Refs. 5, 11, and 12. 

2.4 Discussion 

The present analysis of the Reynolds stress transport equation, based on the 
Yakhot-Orszag renormalization group and (tensorial) 77-expansion summation as 
suggested by Yakhot et al. 9 , has led to a model transport equation incorporating the 
well-known LRR and Rotta models. The analysis gives theoretical support both to 
these models and to the constants sometimes used with them. More significantly, it 
exhibits the LRR and Rotta models as lowest order approximations, and therefore 
also supports their replacement with higher order nonlinear models which would 
be deduced by more accurate approximate summations. The consistency of the 
analysis with higher order effects like the unequal relaxation rates of shear and 
normal stresses has been discussed. 


3. Future Plans 

The nonlinear eddy viscosity representation of the Reynolds stresses 
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T~ij — UiUj ^ K6{j 


c^r 


dUj 

H J 

\dxj dx 


f) 




dUj dUj 
dx p dx p 




dUi dUp 
dx p dxj 


dUj dU p \ (°) f dU v dU p 

dxp dxi ) t3 V dxi dxj 



in which K is the turbulence kinetic energy, e is the dissipation rate, C u , C T \ , C T 2 , C T 3 
are constants, and (0) denotes deviatoric part, was proposed by Yoshizawa 11 in or- 
der to model normal stress effects in shear flows by means of an explicit formula 
for the stresses. The significance of this formula is not limited to this property: 
Yoshizawa’s derivation using a special perturbation expansion, the two-scale direct 
interaction approximation, showed that the expansion could be continued to any 
order in the mean velocity gradient and thereby exhibited the Reynolds stress ten- 
sor as the result of an infinite number of increasingly complex interactions between 
the mean velocity field and turbulence. Related expansions are given in Refs. 5, 12. 

The infinite expansion which contains Eq. (15) can be written symbolically as 


r = KAqS 0 + ^—A\S l -I — y 53 AaSj -I + 

£ £ i<N 2 


R n+1 

e n 


£ AinS? + • • • (16) 

i<N n 


Eq. (16) can be considered a decomposition of the Reynolds stress 

r = r° + r 1 H 


(17) 


where 


r" = £ A in S? < 184 ) 

i<N n i<N n 

where 5 ” denotes a symmetric tensor homogeneous of degree n in the mean velocity 
gradient VU and its transpose, N n denotes the number of linearly independent 
terms of order n (so in Eq. (15), N 2 = 3), and the Ai n are constants. 

In an analysis of the Reynolds stress transport equation by renormalization group 
techniques 35 , we found analogous expansions for the term which governs the return 
to isotropy, 




rK- 




i<N 2 


and for the rapid pressure-strain term 


(19) 


n = -K s 1 + 

5 £ 


s 1 vc/ r + ws 1 ] 
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+ D X 
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+1? H [ Ci2 VC/T + Vt/5?) (0) + Di 2 (SfVU + VU T Sf) 


+ 


( 20 ) 


i<N 2 


Suppose that these sums are introduced into the Reynolds stress transport equation 

f = n' + n- P + £> (21) 


where the dot denotes convective derivative, P = tVU t + VU r is production and D 
denotes the diffusion term. In order to obtain a Reynolds stress transport equation, 
it is necessary to express the sums (19) and (20) in terms of r and VU. Although 
the coefficients A, B, C, D can be explicitly exhibited to any order in perturbation 
theory, they do not have an obvious law of formation. Therefore, the sums (19) 
and (20) can only be approximated by polynomials in r and VU if some hypotheses 
relating the coefficients is introduced. There is no unique hypothesis of this sort, 
but the simplest 35 seems to be 

Bin/Bi = Ai n /A\, n > 2 (22) 

Cin/Ci = D in jD\ = Ai n /Ai, n> 2 

which leads to the Rotta return to isotropy model and to an LRR model for the 
rapid term. The approximation expressed by Eq. (22) can be compared to the 
summation introduced in an analogous context by Yakhot et al. 9 , and to Pade 
approximation: it agrees with the perturbation theory of Eqs. (19), (20) to lowest 
order, but includes effects of all order in VU. 

By evaluating more terms of the perturbation series explicitly and introducing 
an equation like (8) for coefficients of higher order, a hierarchy of models could be 
generated. They would initially be nonlinear in VU, as advocated by Speziale 33 , 
but one might introduce the perturbation series for r • r to obtain a model nonlinear 
in T 1,2 . However, the close analogy between Eqs. (19) and (20) and Yoshizawa’s 
expansion (16) suggests a different approach: namely, use Eqs. (17) and (18) to 
replace 5" in Eqs. (19) and (20) by r"/Aj n . Substitute these modified expressions 
and Eq. (17) into the transport equation Eq. (21), treat r/ 1 as having order | VU | n , 
and separate the terms of like order in | VU | in the standard perturbation theoretic 
fashion. The result is that the terms t" in the decomposition (17), (18) themselves 
satisfy coupled transport equations. 

For simplicity, let us write an approximate system for the r n instead of for the 
r” and assume the most elementary scalar diffusion model. Then the the single 
transport equation for r would be replaced by a system 

f n = -C£r n + Cf (r 

JyT 2 

+C 2 n (r n_1 Vt/ + Vf7 T r n-1 ) (0) + C£ — V 2 T n ,n > 1 (23) 

Since r = , the system (6) should be constrained to satisfy Crow’s condition 

and to contain the exact production term following summation over n. Making the 


n ~ 1 VU T + Vt/r n_1 ) 
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coefficients C R , C", C%, Cp, independent of n reduces the system to a model of the 
LRR form for r = ^ r n . 

This conclusion also follows from Leslie’s analysis 34 of the direct interaction (DIA) 
equations for shear flow. Leslie suggested a perturbative solution for the (tensor) 
correlation function and Green’s function 


u = u° + u 1 + 

G = G° + G 1 + 


(24) 


where U n and G n are of the order | Vt/ | n , and observed that this expansion is 
simultaneously an expansion in powers of the mean strain, and a decomposition into 
symmetry types of increasing complexity. This is also a feature of the expansion 
(16). Substitution of Eq. (20) into the equations of the direct interaction approx- 
imation gives a coupled system for the U n and G n in standard fashion. Then in 
principle, by integrating each equation of this system over all wavenumbers and 
introducing the definitions 


r n = J U n d k 

we could attempt to obtain coupled transport equations for the r n . Unfortunately, 
the derivation of equations for single-point quantities from DIA is not entirely 
straightforward, and more heuristic methods like two-scale DIA 8,14 and renormal- 
ization group are required. 

The simplest model of this type is a two component model, 

t-^KI = t x +t 2 

with 


f 1 =- C 1 r t 1 — —K (VU + VC/ T ) + (C{ - 1) (r 2 W T + VUr 2 ) (0) 

15 

-f- C*2 (t 2 V£/ + WU t t 2 ) (0) + C 1 d V~Vt 1 
f 2 = - C 2 r t 2 + (Cl - 1) (t x VU t + Vt/r 1 ) ( ° ) 

+ C\ (r 1 VU + VU t t x ) (0) + C 2 d V ^ Vt 2 


(25) 


The appearance of r 2 in the equation for r 1 is required by the production term 
constraint which is violated by direct truncation of the system (20) at the second 
order. The properties of this model for simple shear flow, in which dUi/dxj = 
S6*6 j2 follow by setting 


r 1 = ri 2 
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Til 

1 

T 2 = 

T22 
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T 33 _ 


(26) 
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Then C R describes the return to isotropy of the shear stress, and C R the return 
to isotropy of the normal stresses. These relaxation rates can be unequal in this 
theory, an effect predicted by Weinstock’s analysis 3 of the spectral dynamics of the 
return to isotropy. Weinstock suggested further that the individual normal stresses 
relax at different rates: this effect is not accommodated by the present model, but 
could occur in a model in which r 2 is divided into the three tensor components 
t 2 ,t 2 2 ,t 2 of Eq. (18). 

The inequality of the Rotta constants in the present model can be used to over- 
come a defect of models of the LRR type, that in semi- homogeneous flows in which 
a,ij = /K is approximately constant, the ratio an is too small whereas is 
too big 10 . Following Speziale 19 , write Eq. (21) as a system of equations for the 
ratios a^- and set = 0. There results, 


2ai2P — ( C l R — l) — 


(C R — l) — 4p 2 (C 1 — l) 022 + C 2 a n 



an 


a 22 = 


( -C 2 r + l) + va 12 

[l (<?i ~ !) ~ t^]rm2 
(~C R + 1) + 77012 


When P/e = —77012 = 1, 


a 12 — 


4_J_ 
15 C R 


C\-l 

Cr 


Cl 

022 — TTf On 


(27) 


Oil — 

o 22 = 


-f (cl - 1) + \cl\icl 
I (Cl - 1) - \cl\cl 


(28) 


Eq. (28) shows that setting C R < C R both increases an and decreases ai 2 and 
thus improves the agreement between theory and experiment. The trend required 
here is consistent with Weinstock’s findings 3 that the Rotta constant for the normal 
stresses should be smaller than the shear constant and should take values close to 
1 . 0 . 

The inequality of the coefficients describing the rapid terms affects the behavior 
of the model under rapid distortion. The rapid distortion analysis of passively 
strained turbulence predicts that a one-component limit state is reached in which 
nj = 2K8ji6ji, 77 = SK/e — > 00 and P/e is finite 21 . Whether or not this limit 
can in fact occur as an asymptotic state, a stress model should accommodate it 
because it can exist approximately in steady flows as a ’’spatial transient” in strongly 
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inhomogeneous regions of very high convective or diffusive transport 17 . Moreover, 
models of the LRR form do not capture the transient evolution of rapidly distorted 
flows well. Let us assume that all model coefficients are functions of the ’’state” 
of turbulence, following Shih and Lumley 2 ; the precise parametrization of the state 
will be left to future investigations. From Eq. (13), it is evident that P/e = —77012 
will be of order 77 in the one-component state in which an = 4/3,022 = —2/3 unless 


-I ~ ') + H 


4_ 

15 


in this state. But Eq. (13) also shows that On = 4/3,022 = —2/3 requires 


Cl ~ 0 , Cl ~ 0 

in this state. These conditions are inconsistent in a model in which C\ = Cl and 
C\ = Cl, but are clearly consistent with the present proposal. 

The advantages of this model are consistency with a systematic perturbation 
theory, the possibility of unequal relaxation rates for normal and shear stresses 
in relaxing strained turbulence, improved agreement with experimental data for 
universal ratios in simple shear flows, and the possibility of accommodating the 
one-component limit of rapid distortion theory. 
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Bypass Transition in 
Compressible Boundary Layers 

J. J. Van der Vegt 
1. Motivation and Objectives 

Transition to turbulence in aerospace applications usually occurs in a strongly dis- 
turbed environment. For instance, the effects of free-stream turbulence, roughness 
and obstacles in the boundary layer strongly influence transition. Proper under- 
standing of the mechanisms leading to transition is crucial in the design of aircraft 
wings and gas turbine blades, because lift, drag and heat transfer strongly depend 
on the state of the boundary layer, laminar or turbulent. Unfortunately, most of 
the transition research, both theoretical and experimental, has focused on natural 
transition. Many practical flows, however, defy any theoretical analysis and are ex- 
tremely difficult to measure. Morkovin 5 introduced in his review paper the concept 
of bypass transition as those forms of transition which bypass the known mecha- 
nisms of linear and non-linear transition theories and are currently not understood 
by experiments. 

In an effort to better understand the mechanisms leading to transition in an 
disturbed environment, experiments are conducted studying simpler cases, viz. the 
effects of free-stream turbulence on transition on a flat plate, Sohn and Reshotko 14 
and Wang et al. 19 . It turns out that these experiments are very difficult to conduct, 
because the generation of free-stream turbulence with sufficiently high fluctuation 
levels and reasonable homogeneity is non trivial. For a discussion see Morkovin 5 . 
Serious problems also appear due to the fact that at high Reynolds numbers the 
boundary layers are very thin, especially in the nose region of the plate where the 
transition occurs, which makes the use of very small probes necessary. 

The effects of free-stream turbulence on transition are the subject of this re- 
search and are especially important in a gas turbine environment, where turbu- 
lence intensities are measured between 5 and 20%, Wang et al. 19 . Due to the fact 
that the Reynolds number for turbine blades is considerably lower than for aircraft 
wings, generally a larger portion of the blade will be in a laminar-transitional state. 
Turner 15 shows that the effect of free-stream turbulence on transition significantly 
increases when the free-stream turbulence levels become larger than 5% and is ac- 
companied with a significant increase in heat transfer. Recently Rai and Moin 11 
presented a direct numerical simulation of transition to turbulence on a flat plate 
in a free-stream with Mach number .1 and turbulence levels at the leading edge 
of about 2.75%. Direct numerical simulations offer a unique opportunity to study 
specific phenomena, while excluding disturbances from other sources. The com- 
putations from Rai and Moin show some impressive results, especially regarding 
intermittency and turbulent spots. Their numerical simulation, however, has the 
same problem as with most of the experiments, namely a very low Mach number, 
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while many applications operate in the transonic regime. Due the nature of their 
numerical scheme, a non-conservation formulation of the Navier-Stokes equations, 
it is a non-trivial extension to compute flow fields in the transonic regime. 

This project aims at better understanding the effects of large free-stream tur- 
bulence in compressible boundary layers at Mach numbers both in the subsonic 
and transonic regime using direct numerical simulations. The present project aims 
at computing the flow over a flat plate and curved surface. This research will 
provide data which can be used to clarify mechanisms leading to transition in an 
environment with high free stream turbulence. This information is useful for the 
development of turbulence models, which are of great importance for CFD appli- 
cations, and are currently unreliable for more complex flows, such as transitional 

flows. 

2. Accomplishments 

Direct simulations of transition in compressible flows with both shocks and bound- 
ary layers requires an extremely accurate and efficient scheme. Several conflicting 
requirements present a serious challenge which cannot be met by existing numerical 


• The small grid spacing in the boundary layer makes an implicit scheme necessary, 
because an explicit scheme would have a severe time step limitation. Implicit 
schemes usually are not time accurate and rather dissipative. 

• Higher order accurate schemes are necessary but higher order accurate schemes 
generally do not give non-oscillatory solutions around discontinuities, such as 
shocks. Many of the popular non-oscillatory shock capturing schemes, such as 
TVD (Total Variation Diminishing) methods, are only first order accurate in 
multi-dimensional flows and even in one-dimension they reduce to first order at 
non-sonic local extrema. 

In order to satisfy these conflicting requirements a significant effort has been made 
to improve and combine several successful numerical schemes. A fully implicit and 
time accurate code for the solution of the three-dimensional compressible Navier- 
Stokes equations in general geometries has been written and tested. Higher order 
accuracy and shock capturing are implemented using an Essentially Non-Oscillatory 
(ENO) scheme. Time accuracy is obtained using a Newton method. 

In the next section a brief description of the numerical scheme will be given 
followed by the discussion of a series of tests aimed at validating the code. 

2.1 Numerical Scheme 

The compressible Navier-Stokes equations are solved using a finite volume method .| 
A detailed discussion of finite volume and difference methods can be found in 
Vinokur 18 . The integral formulation of the Navier-Stokes equations, assuming all 
variables are continuous in time, is given by: 


schemes: 



( 2 . 1 . 1 ) 
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Here V(t ) and surface S(t) are the volume and outer surface of the domain Q 
and n an outward unit normal vector at S. The vector U represents the conserved 
variables: ( p , pu, pv, pw , e) T , with p the density, u, v and w the velocity components, 
and e the total energy. The tensor P is defined as P = S + V, with 5 the inviscid 
contribution defined as: 


/ pu \ 

/ pu \ 
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and V the viscous contribution: 


V i = 



V 2 = 



V, = 


The shear stress tensor T, with components ( 71 , 72 , 73 ) is defined as: 


( 2 . 1 . 2 ) 



(2.1.3) 


T = -J-(/z(Vu + Vu T ) + A(V • u)2) 
Re 


(2.1.4) 


and the heat flux q as: 

/n 1 r\ 

q " ( j - l )ReM^Pr 

The variables p, T, p, A and k represent the pressure, temperature, first and 
second viscosity coefficient and thermal conductivity, respectively. The coefficients 
Re, M^, and Pr are the Reynolds, Mach, and Prandtl numbers. All variables are 
non-dimensionalized using free-stream variables and a characteristic length L. 

The Navier-Stokes equations are solved using a finite volume method because we 
seek a weak solution in order to capture shocks in high Reynolds number flows. 
The finite volume method is also the most natural way to satisfy the conservation 
properties of the differential equations. After subdividing the volume V into a set of 
disjunct cells we obtain the finite volume discretization for a cell with index i, j, k: 


~ (vu«,») + F! + j ,i, k - t'l-i.i.k + F? J+ i,k - 
ol ( 2 . 1 . 6 ) 
where a barred quantity with index i,j,k is an average of the unbarred quantity 
over the cell with index i,j, k and indices i ± j ± % and k ± | refer to values at 
the cell faces. The fluxes F* at the cell faces are defined as: 


+ F? 




P = S l -F 


(2.1.7) 


with S l the cell face in the direction i. The computation of the cell face S l and 
volume V has to be done with great care in order satisfy the geometric conservation 
law, for details see Vinokur 19 : 
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Flux Approximation 

The crucial part in the development of a finite volume method is the approxi- 
mation of the fluxes at the cell faces. The flux Fj +i is computed using the Osher 
approximate Riemann solver. The first order accurate conservative flux is given by: 


F 1 

r t+ 




| du) 


( 2 . 1 . 8 ) 


with equivalent expressions for the other two directions. The integral is computed 
along a path in phase space, connecting the points with index (i,j, k ) and (i+lj, k). 
Along each subpath a Riemann problem is solved, which is used to determine the 
intermediate states. In this way exact expressions can be derived for the path inte- 
grals. More details about the implementation of the Osher scheme can be found in 
Osher and Solomon 6 , Osher and Chakravarthy 7 , Chakravarthy and Osher 1 and Rai 
and Chakravarthy 10 . The Osher approximate Riemann solver is the most accurate 
approximate Riemann solver and satisfies the entropy condition, contrary to the 
Roe approximate Riemann solver which needs an entropy fix to eliminate steady 
expansion shocks. The Osher scheme captures steady shocks in at most two points. 

The most important reason for the choice of the Osher scheme, however, has 
been its low numerical dissipation in boundary layers, Koren 4 . Most schemes for 
the Euler equations are very dissipative in the boundary layer and not well suited 
for direct numerical simulations. In earlier work, Van der Vegt 17 , modifications to 
flux vector splitting schemes were discussed to alleviate this problem, but although 
significant improvement was achieved on steady laminar boundary layers, it was 
not possible to reach accuracy levels necessary for direct simulations. 

Higher order spatial accuracy 

Direct simulations require a high accuracy which cannot be achieved with stan- 
dard second order schemes. It is fairly straightforward to derive higher order accu- 
rate finite difference schemes, but shock capturing then will not be possible. The 
development of higher order accurate, multi-dimensional finite volume schemes, 
capable of shock capturing still is an area of active research, and has been an im- 
portant subject in this project. A significant effort has been made to combine the 
Osher approximate Riemann solver, discussed in the previous section, with an ENO 
scheme. Results of this work are described in Van der Vegt 17 , where the different 
ENO methods are discussed and results of various tests are discussed. 

Higher order accuracy of a finite volume method can be defined is various ways. 
One approach is to define higher order accuracy with respect to the cell averaged 
values. This resembles most closely the finite volume description, which gives equa- 
tions for the cell averaged values. Another definition of higher order accuracy uses 
the point values at the cell center, which is used in conservative higher order finite 
difference methods. Both approaches are being used. For subsonic flows currently 
the fifth order scheme, developed by Rai 9 , is used, which is based on a Taylor series 
expansion of the flux vector along the lines presented by Osher and Chakravarthy 8 . 
This method is a conservative finite difference scheme. It has the benefit that it is 
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simple to implement in more dimensions, but does not allow shock capturing. Os- 
her and Chakravarthy 8 demonstrated how to make these schemes TVD and allow 
shock capturing, but they are not very useful and only first order accurate globally. 
The scheme therefore is used in its unlimited form, limiting its application to flows 
without discontinuities such as shocks. The scheme also is rather expensive and 
work is in progress to improve its efficiency. 

For flows with shocks research has been carried out to develop higher order ac- 
curate ENO schemes. ENO schemes use an adaptive stencil where a searching 
algorithm tries to find that part of the flow field surrounding a cell which is the 
smoothest. Then a conservative, higher order accurate interpolation method is used 
to ’’reconstruct” the point values from the cell averaged values. Due to the fact that 
the interpolation process only uses data from the smooth part of the flow field nu- 
merical oscillations will be minimized. In this way uniform higher order accuracy 
can be obtained. The first ENO methods were developed by Harten et al. 3 , and 
later modifications were proposed by Shu and Osher 12,13 . In Van der Vegt 17 the 
different methods were compared and it was found that the ENO scheme, using 
primitive function reconstruction from cell averaged variables with the Cauchy- 
Kowalewski procedure for time integration combined with the Osher approximate 
Riemann solver, was the most accurate and robust. In one dimension it has been 
successfully used up to fifth order accuracy, but due to the fact that in multi- 
dimensional flows currently dimension splitting is used, its accuracy is limited to 
second order in more than one dimension. Work to extend this scheme to higher 
order accuracy in multi-dimensional flows is in progress. 

Time integration 

Due to the very small gridspacing necessary at the wall and in critical layers 
explicit time integration would result in a serious time step limitation. To alleviate 
this problem implicit time integration has to be used, but most implicit time inte- 
gration schemes make assumptions in the implicit part which reduce or eliminate 
time accuracy. The development of implicit and time accurate numerical schemes 
therefore has been a significant part of this research. Time accuracy is obtained 
using the Newton method discussed in Rai 9 , which solves the non-linear system of 
equations in the implicit time integration scheme using a Newton method. Rai uses 
this method also to reduce the error caused by approximate factorization. We do 
not use approximate factorization but solve the whole matrix system iteratively, see 
Van der Vegt 16 , and need the Newton scheme only to reduce the error due to the 
time linearization. This iterative scheme also has the benefit that it is not neces- 
sary to use an exact linearization of the flux vectors, which can be very difficult 
and time consuming to obtain. First order Steger- Warming flux vector splitting is 
used in the implicit scheme, while a higher order accurate spatial discretization is 
used for the explicit part. At each time step the Newton iteration is performed 
such that the accuracy of the higher order explicit part is maintained. The use of 
an approximate linearization, however, limits the maximum time step and work is 
in progress to evaluate if a more accurate linearization would improve the perfor- 
mance and robustness of the scheme. Especially at high Mach numbers there still 
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are convergence problems for large Courant numbers when the scheme is used to 
obtain steady state solutions. 

2.2 Results and Discussion 

Several computations have been carried out to test the code and the various 
numerical schemes. In order to test the Osher approximate Riemann solver and 
the ENO schemes a large number of shock tube calculations have been carried out, 
a detailed description of this work can be found in Van der Vegt 17 . The different 
ENO schemes tested were the ENO method from Harten et al. 3 , using primitive 
function reconstruction and either Runge-Kutta time integration or the Cauchy- 
Kowalewski procedure, and the Shu and Osher flux- ENO scheme. In all cases the 
Osher approximate Riemann solver was used and the effect of the ordering of the 
eigenvalues, viz. natural and reversed ordering, has been investigated. Four cases 
with different difficulties were tested, see Table 1. The performance of the different 
schemes was reasonable in most cases, but it turned out that the ENO scheme 
with primitive function reconstruction and the Cauchy-Kowalewski procedure for 
time integration (ENO-CK) was the most accurate and robust. Some of the results 
obtained with this method, are shown in Figures 1 and 2. Figure 1 shows a left 
moving shock followed by a contact discontinuity and a right moving expansion 
wave. A difficult problem for the ENO schemes in case A is the fact that in the 
initial stages there are not enough grid points available in the region between the 
shock and contact or shock-shock. The ENO scheme searches for the stencil which 
gives the smoothest part of the flow field around a grid point or cell. In these 
cases there exist in both directions a discontinuity and there are not enough points 
available to build a non-oscillatory higher order reconstruction. This problem exist 
for all ENO schemes but the ENO Cauchy Kowalewski scheme is the least sensitive 
for it. The other methods have mild to strong oscillations in these areas. One way 
to solve this problem is to reduce the accuracy locally till enough grid points are 
available to create a higher order reconstruction, but this is a problem which still 
needs further attention. Case B, Figure 3, shows a left moving expansion wave and 
a right moving contact discontinuity and shock. One of the problems in this case is 
the appearance of a sonic point, which gives a small jump of 0( Ax) at first order. 
The shock tube tests showed that it is possible to use a higher order scheme for 
flows with discontinuities, but the convergence of these higher order schemes can 


Case 

Pl 

Pr 

U L 

Ur 

T l 

Tr 

A 

15000 

98400 

0 

0 

1378 

4390 

B 

988000 

9930 

0 

0 

2438 

2452 

C 

10000 

100000 

0 

0 

2627 

272 

E 

573 

22300 

2200 

0 

199 

546 


Table 1. Initial Conditions Shock Tube Tests. 
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Fig. 1. Case A, density at t = .2, 2nd and 5th order ENO-CK. 




Fig. 2. Case B, density at t = .2, 2nd and 5th order ENO-CK. 


be at most first order around these discontinuities. Tests of all the ENO schemes 
on smooth solutions showed that they all reached the proper level of accuracy. Test 
are currently underway to check if the ENO schemes give higher order accuracy in 
regions outside discontinuities as they are supposed to. This is an important test 
to see if these methods are capable of shock-turbulence interaction simulations. 

In order to test the shock capturing properties of the code the flow field around 
a circular cylinder at Mach 8 has been computed. Although the flow field was two- 
dimensional the three-dimensional code was used to check if the flow field remained 
exactly two-dimensional and the geometric conservation law was satisfied. Figure 
3. presents contours of the pressure at a spanwise station along the cylinder. The 
solution is the same at all stations. The flow field consists of a strong bow shock 
where the Mach number changes from 8 to about 2.8 behind the shock at the 
symmetry line. Apart from the strong shock another aspect of this case is the 
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Figure 3. Pressure field of flow around a cylinder at Mach 8 


fact that the flow field in the stagnation region ahead of the cylinder is subsonic. 
The sonic line is at about 45 degrees with the flow angle and a smooth transition 
is observed from the subsonic to the supersonic region. This case has also been 
computed by Osher and Chakravarthy 7 and the results compare well. To test the 
ability of the code to simulate transitional flows which is a crucial test before bypass 
transition can be simulated currently computations are done on natural transition 
in a flat plate boundary layer. A comparion is made with the results of Fasel et 
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Figure 4. Vertical velocity field in a flat plate boundary layer at Mach 0.08 with 
periodic suction and blowing (vertical direction enlarged 20 times). 


al. 2 , which computed transition in an incompressible boundary layer. In order to 
make the comparison as accurate as possible a very low Mach number, & > = .08 
was chosen. This Mach number is approximately the lower limit for the numerical 
scheme and despite the fact that the computations are fully implicit a severe time 
step limitation is imposed by the sound waves. To start the simulation first a steady 
laminar boundary layer is computed, which is a non-trivial task because a very high 
accuracy is needed in this computation. The disturbances at the beginning of the 
plate have an amplitude of 10 -4 and transition simulations require a very low nu- 
merical dissipation. The disturbances are generated by periodic suction and blowing 
in a strip somewhat downstream of the inflow boundary. This is done because there 
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are always numerical disturbances due to the fact that the inflow boundary condi- 
tions are not perfect. The disturbances are generated in a region of the boundary 
layer which is linearly stable. When they move downstream first the transients are 
damped and after they move in the unstable region the Tollmen Schlichting waves, 
which are the most unstable two-dimensional waves, are amplified. In order to ac- 
commodate for the fact that the subsonic outflow boundary conditions, which are 
essentially inviscid, are not perfect in a boundary layer a buffer layer is added to the 
plate to damp as much as possible the reflections coming from this boundary layer. 
This is the same procedure as used by Rai and Moin 11 . The number of grid points 
in this computation is 340 x 82. It turned out that the most efficient way to obtain 
the steady boundary layer was to first start with the second order scheme running 
at a Courant number of 120 using the implicit Euler time integration and to switch 
to the fifth order scheme after most of the transients are damped out. The fifth 
order scheme has a very low dissipation and would otherwise take a long time to 
converge. The maximum CFL number for which the fifth order scheme is stable is 
approximately 160. Ater the steady state was obtained periodic suction and blowing 
were added and the fifth order scheme was used with the Newton time integration 
scheme at a CFL of 60. Figure 4. shows a preliminary result of this computation 
and it clearly shows the gradual build up of the boundary layer instability. The 
results are currently analysed to make a comparison with the incompressible results 
of Fasel et al. 2 . 

3. Future Plans 

Further testing of the code will have to be done for more complicated cases. 
Currently a comparison with the results of stability of a flat plate boundary layer 
at low Mach number with the results from Fasel et al. 2 is being completed. If 
this comparison and equivalent tests for high Mach number boundary layers are 
satisfactory a simulation of bypass transition on a subsonic flat plate boundary 
layer will be made, followed by simulations of a boundary layer on a curved plate. 
Also work has to be continued to develop higher order accurate ENO schemes for 
multi-dimensional flows. This is crucial for direct simulations of transonic flows 
with both shocks and turbulence. 
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Runge-Kutta Methods Combined 
with Compact Difference Schemes 
for the Unsteady Euler Equations 

Sheng-Tao Yu 


1. Motivation and Objective 

Recent development using compact difference schemes to solve the Navier-Stokes 
equations show spectral-like accuracy 1 ’ 2 . In this paper, we report further study 
of the numerical characteristics of various combinations of the Runge-Kutta (RK) 
methods and compact difference schemes to calculate the unsteady Euler equations. 
Conventionally, the accuracy of finite difference schemes is assessed based on the 
evaluations of dissipative error. The objectives are reducing the numerical damp- 
ing and, at the same time, preserving numerical stability. While this approach has 
tremendous success solving steady flows, numerical characteristics of unsteady cal- 
culations remain largely unclear. For unsteady flows, in addition to the dissipative 
errors, phase velocity and harmonic content of the numerical results are of concern. 
As a result of the discretization procedure, the simulated unsteady flow motions 
actually propagate in a dispersive numerical medium. Consequently, the dispersion 
characteristics of the numerical schemes which relate the phase velocity and wave 
number may greatly impact the numerical accuracy. The objective of the present 
paper is to assess the numerical accuracy of the simulated results. To this end, 
the Fourier analysis is performed to provide the dispersive correlations of various 
numerical schemes. 

First, a detailed investigation of the existing RK methods is carried out. A 
generalized form of an N-step RK method is derived. With this generalized form, 
the criteria are derived for the three and four-step RK methods to be third and 
fourth-order time accurate for the non-linear equations, e.g., flow equations. These 
criteria are then applied to commonly used RK methods such as Jameson’s 3- 
step and 4-step 3 ’ 4 schemes and Wray’s algorithm 5 to identify the accuracy of the 
methods. For the spatial discretization, compact difference schemes are presented. 
The schemes are formulated in the operator-type 6 to render themselves suitable for 
the Fourier analyses. The results of the analyses provide CFL limits, the numerical 
dispersion relations, and the artificial damping required for stable and time-accurate 
solutions. 

Finally, the performance of the numerical methods is demonstrated by numeri- 
cal examples. The first case is a quasi-one-dimensional calculation of the acoustic 
admittance in a converging nozzle. The CFD results are compared with Tsien’s 
analytical solution 7 ; the harmonic content of this flow field is limited to one fre- 
quency mode. All numerical schemes of concern provide accurate solutions. The 
second case is a one-dimensional simulation of a shocked sound wave. The harmonic 
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content is complex and distinct differences between various schemes are observed. 
The results are also compared with the analytical solution provided by Morse and 
Ingard 8 . In the one-dimensional cases, details of the numerical methods in setting 
up the initial conditions and the perturbation on the computational boundary are 
described. 

The third case is a two-dimensional simulation of a Lamb vortex 9 in an uniform 
flow. This calculation provides a realistic assessment of various finite difference 
schemes in terms of the conservation of the vortex strength and the harmonic content 
after travelling a substantial distance. The numerical implementation of Giles’ non- 
reflective equations 17 coupled with the characteristic equations as the boundary 
condition is discussed in detail. Finally, the single vortex calculation is extended 
to simulate vortex pairing 10 . For the distance between two vortices less than a 
threshold value, numerical results show crisp resolution of the vortex merging. 

2. Work Accomplished 
2.1 Numerical Method 

The Euler equations in Cartesian coordinates can be cast into a vector form: 


where Q is the unknown vector and E; is the inviscid flux in the Xi direction. 
The Runge-Kutta algorithm is applied as the temporal discretization and the sec- 
ond, fourth, and sixth-order compact difference schemes are applied to the spatial 
discretization. 

2.2.1 The Runge-Kutta Method 

The use of the Runge-Kutta methods for flow equations stems from the applica- 
tion of the methods to solve ordinary differential equations (ODEs). An ODE has 
one independent variable and its solution is obtained by integrating the equation 
from its initial condition. When one applies the Runge-Kutta method to the flow 
equations, time is treated as the independent variable as is in an ODE, and the 
convective terms are taken as the inhomogeneous part of the equations, such as 


Notice that the boldface symbols which represent vectors have been temporarily 
dropped for typographic convenience. In addition, all results in the following dis- 
cussion are valid for both scalar and vector equations. 

The Runge-Kutta methods have algorithms of the form 



( 1 ) 



( 2 ) 


QU+l = QU _J_ /\tR(Q n ,At), 


( 3 ) 
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where the increment function R(Q n ,At) is a suitable chosen approximation to the 
inhomogeneous part of the equation, that is, R(Q). In general, the calculation of 
the increment function R is subdivided into N steps on the interval t n < t < t n+1 . 
And the final increment function R is a weighted average of the inhomogeneous 
terms evaluated at the different steps on the interval t n <t < t n+1 , that is 

Q l = Q n + At(a n R n ), 

Q 2 = Q n + At(a 2 l R n + a 22 R l ), 

Q 3 = Q n + At(a 31 R n + a 32 R 1 + OC 33 R 2 ), (4) 


Q n+1 = Q n + At(a m R n + a^R 1 + • • • + a NN R N ~ 1 ), 

where the superscript n, 1 , 2 ,..., and n+ 1 denote the time steps on the time interval 
t n < ti < t 2 < ... < tjv < t n+1 , and otij is the weighting factor for the step i and 
term j. There are Y^iLi * weighting coefficients to be determined and an infinite 
number of coefficient sets can be chosen. However, certain criteria must be met for 
the algorithm to retain high-order accuracy. 

In what follows, the criteria of the coefficient set of a 3-step Runge-Kutta method 
to be third-order accurate is given. To proceed, we follow the conventional approach 
and expand all inhomogeneous terms R l in Eqn. (4) to a Taylor’s series about R n 
and drop all terms in which the exponent of At is greater than 3. The result 
is compared with its analytical counterpart by equating terms in like powers of 
At. The result is tabulated in Table 1. For the convenience of the discussion, the 
following simplification of symbols is activated: R denotes R(Q n ), Q denotes Q n , 
R' denotes ( dR/dQ) n , and R" denotes ( d 2 R/dQ 2 ) n . In addition to the equality of 
the coefficients of all the powers of At, we also want the equality of the coefficients 
of the functions of R, R' , and R" . As a result, we find the criteria of the coefficients 
for the 3-step RK methods to be third-order accurate as, 

#31 T 0:32 + <*33 = 1 ) ( 5 fl) 

Qui<*32 + <*33 (<*21 + <* 22 ) = 2> (5h) 

anC*32 ( a21 + q 22) 2q: 33 = (5c) 

<*110=220:33 = g- ( 5 d) 

Equations (5a) and (5b) are the criteria of first and second order accuracy, re- 
spectively. The remaining equations are of third order term. Since four equations 
contain six unknowns, the system is underdetermined, and two of the coefficients 
may be chosen arbitrarily. The obvious choice is to let the two coefficients be null 
to reduce intermediate storage and numerical operations. According to Eqn. (5d), 
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none of an where i = 1,2,3, could be zero, and one can set the two of the three re- 
maining coefficients to be zero. Therefore, at least one of the intermediate steps has 
two non-zero coefficients. Consequently, one needs to store two steps of intermediate 
solutions for the 3-step, third-order RK methods. 

A 3-step Runge-Kutta method proposed by Jameson et al. 3 to solve flow equations 
is 


Q 1 = Q n + A tR n , 

Q 2 = Q n + y( r " + r ">’ ( 6 ) 

It can be shown that the weighting coefficients of the 3-step method satisfy only 
Eqns. (5a) and (5b). And the method is second-order accurate in time. 

Wray 5 proposed another 3-step method, 

Q 1 = Q n + At(^iT), 

15 

Q 2 = Q 1 + At(^R" - ^R 1 ), (?) 

Q n+1 =Q 2 + At(^R n - ^H 2 ). 

This formula may be manipulated to fit the generalized form as proposed in Eqn. 
(4), and we obtain, 

Q 1 - Q n + At(-^R n ), 

15 

Qi =Qn + At(±R n + ^R 1 ), (8) 

Q n +1 = qn + Af ( 1 R n + ^R 2 ). 

Wray’s coefficients match all the equations in Eqn. (5) and therefore the scheme is 
third-order accurate. In this formula, a 32 is set to zero and two sets of solutions 
are needed in the second and the third steps. The calculation can be carried out by 
either the vectorized algorithm proposed by Wray, or straightforward calculation 
according to Eqns. (7) and (8). 

A similar procedure can be applied to the 4-step Runge-Kutta methods, and the 
criteria for the scheme to be fourth-order accurate are: 


&41 + &42 + C*43 + #44 — 1, 
OUl £*42 + (<221 + <222)<243 + (<231 + O32 + <233)0:44 = -, 
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0^42 + (<*21 + <*22) 2 <*43 + (<*31 + <*32 + <*33) 2 <*44 = (® c ) 

<*11 <*22 <*43 + [<*11 <*32 + (<*21 + <*22 )<*33] <*44 = g j (^<0 

<*11 <*42 + (<*21 + <*22 ) 3 <*43 + (<*31 + <*32 + <*33) 3 <*44 = T> (^ e ) 

Ct22 <*43 + (<*21 + <*22 )<*11 <*22<*43 
z 

+ ~ [(<*21 + <*22 ) 2 <*33 + <*ll <*32] <*44 
+ ( 0:31 + <*32 + <*33) [<*11 <*32 + (<*21 + <*22 )<*33] <*44 = g? (9/) 

<*11 <*22 <*33 <*44 = 24* (9<?) 

Equations (9a) and (9b) are for first and second-order accuracy, respectively. Equa- 
tions (9c) and (9d) are for third-order accuracy. The remaining equations are for 
the fourth-order terms. Here, seven equations contain ten unknowns, and three of 
the coefficients may be chosen arbitrarily. 

A 4-step RK method attributed to Kutta 11 for solving ODEs was adopted by 
Jameson et al. 3 to solve the flow equations. The algorithm can be expressed as, 


Q 1 = Q n + Y Rn ' 

Q* = Qn + ^±R\ 

Q z = Q n + AtR\ 

QU+1 = QU + + 2R i + 2 R 2 + R z ). 


( 10 ) 


The coefficients satisfy Eqn. (9) and the algorithm is fourth-order accurate. How- 
ever, this method requires all four intermediate solutions in the final step. As a 
result, the use of this scheme for large-scale calculations is undesirable. 


Later on, Jameson and Baker 4 proposed another 4-step algorithm, 


Q 1 =Q n + —R n , 

Q 2 = Q n + Y R '' 

Q 3 = Q n + Y R2 ’ 

qtx+ 1 = Qn + 


( 11 ) 


This scheme was proposed to calculate the steady state solutions and the transient 
solutions were not of concern. For that purpose, the algorithm is convenient to 
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program and no intermediate solution needs to be stored. For the present investi- 
gation, however, the weighting coefficients of this method satisfy only part of Eqn. 
(9) and the algorithm is second-order accurate. Nevertheless, this formulation is 
favorable compared to a 2-step, second-order RK method because part of the third 
and fourth-order terms are satisfied, namely, Eqns. (9d) and (9g). Consequently, a 
larger marching step, i.e., a larger CFL number, could be used. 

2.2.2 Compact Difference Schemes 

The remaining task of discretizing the flow equations is the spatial differencing 
of the inviscid fluxes. An effective manner for generating a high-order, central 
difference scheme is the compact difference method. The scheme is obtained by 
using only three and five points to achieve fourth-order and sixth-order accuracy 
in space, respectively. The gain in the accuracy is not based on the involvement of 
more points as in the conventional approach, but on implicitly solving the derivatives 
simultaneously at all locations. According to Hermite’s generalization of a Taylor’s 
series, 12 one can get 


«j_! + 4 u { + u i+1 = ~ Ui ~ i) + 0(Ax 4 ), (12) 

u t _ x + 3 u'i + u i+x = 2 + 28u i+1 - 28 - Ui - 2 ) + 0(Ax 6 ), (13) 

where the superscript ' represents the spatial derivatives. Equation (12) is the 
fourth-order method and Eqn. (13) is the sixth-order one. When the fourth- 
order method is used in the interior nodes, a third-order biased implicit scheme 13 
is adopted for grid nodes at the computational boundary, such as 


2 u x + 4 n 2 = ^ 3 ) 

2 u max + ^ U max-1 ~ ~^(5u m ax ~ ^ u max-l ~ Umax- 2 ) + 0(Ax ). 


(14) 


When the sixth-order method is used in the interior nodes, the fourth-order scheme, 
Eqn. (12), is used at locations one grid node away from the boundary and the 
third-order biased scheme is used at the boundary. The application of the implicit 
compact difference schemes with the appropriate boundary conditions involves the 
inversion of a scalar tridiagonal matrix. The inversion of the matrix incurs little 
penalty in terms of CPU time. 


2.3 Fourier Analysis 

By definition, any function, u(x,t), which is continuous, periodic, and square 
summable can be expressed in a Fourier series expansion at a constant time, 


OO 

u{x,t)= Y, *(M)e i2wfcl/L (15) 

k= — OC 
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where L is the period of the function u(x,t), k is the wave number, and i is \/— I- 
The Fourier coefficient is defined as, 

i rL/% 

u(k , t) = jr / u(x, t)e~ iUkx ' L dx. (16) 

L J-L/2 

In the Fourier analysis of a finite difference scheme, functions are defined at discrete 
points. The discrete Fourier series and its coefficients are defined analogous to their 
continuous counterparts, 


K - 1 

Uj = ^ u n (k)e t2nk ^ K , j = 0, ±1, ±2, • • • , ±oo, 

fc =0 
1 K 

u n (k) = u?e~ i2 * kj/K ’ k = 0, 1, • • • , K - 1. 

3 = 1 


(17) 


Here, the length of the computational domain L is decomposed into K grid nodes 
(L = KAx). The superscript n denotes the time step and the subscript j is the 
spatial location index. Similar to the continuous function, the algebraic system in 
terms of the function expiifl'nkj /K) is periodic over the computational domain L 
(or K) and orthonormal, such as 


e i2*kj/K _ e i2nk(j+K)/K 

1 K ~ l ( 18 ) 

^ e a*k l jlK e -t2*k a} /K = 

3 = 0 

where 6 kl k 2 is the Kronecker delta. Therefore, the establishment of the discrete 
Fourier series and coefficients is self-sufficient, and is not an approximation of its 
continuous counterpart. 

As shown in Eqns. (17) and (18), the harmonic content of the discretized equation 
is limited to the number of grid nodes used in the computational domain. A discrete 
solution uf at a location ( j ) and time (n) is a linear combination of K wave modes. 
The Fourier analysis is performed by substituting each wave mode of the discrete 
Fourier expansion, Eqn. (17), into the discretized flow equations to calculate the 
amplification factor, g(k), which is defined as 


9(k) = 


u n+1 (k ) 
u n (k ) 


(19) 


The procedure is repeated for all wave modes ( k = 0, 1, - - * ,K - 1) and the full 
spectrum of the amplification factor is obtained. In this process, we map the func- 
tion u in terms of spatial variable x on the interval [— L/2,L/2] to the wave number 
space on [- 7 r, 7 r] assuming that the analysis is local for an infinite and periodic do- 
main. Therefore, the solution of the amplification factor on the interval [— 7r, 0] is 
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the complex conjugate of that on [0,7r]. For this reason, the results of our Fourier 
analyses are presented on the interval [0,7r]. 

In the present investigation, one-dimensional equations are considered for the 
Fourier analysis. In addition, we can perform a similarity transformation to trans- 
form the one-dimensional Euler equations to their characteristic form, i.e., three 
decoupled scalar equations. Consequently, a scalar, advective equation on a peri- 
odic domain is adopted as the model equation in our analysis, 




( 20 ) 


where the phase velocity (A) is equivalent to the eigenvalues of the Euler equations, 
namely u - c, u + c, and u where u is velocity and c is the speed of sound. The 
phase speed (A) is treated as a parameter in the Fourier analysis to avoid the 
Fourier convolution, therefore the analysis is linear. For the unsteady calculations, 
the requirement of the time resolution of the flow field restricts the time marching 
step. In other words, the variations of flow properties between time steps are small. 
Thus, linear analysis is a viable tool. 

In what follows, the procedure to obtain the amplification factor of the model 
equation discretized by Runge-Kutta methods and compact differences is illustrated. 
First, the generalized forms of the amplification factor for the 3 and 4-step Runge- 
Kutta methods are derived. These representations of the amplification factors are 
independent of the spatial discretization schemes. From the equations of Wray’s 
3-step scheme, Eqn.(7), we have 


S 2 - 1 + 4 Z+ i2 Zp1 ’ 

9 = l + \ z + l z 9 2 , 

where 

g'Q" = Q\ 

9 2 Q" = Q 2 , ( 22 ) 

gQ “ = 

The variable Z represents the spatial discretization applied to the convective term ( 
—A du/dx ). By substituting the amplification factors of the intermediate steps, g 1 
and g 2 , into the last step of Wray’s algorithm, we obtain the amplification factor, 
g , for the 3-step scheme, 


<j = l + Z + iz 2 + iz 3 . 


(23) 
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It is interesting to note that Eqn. (23) can be directly derived from the Taylor’s 
series expansion by adopting the invariance property of the time and spatial deriva- 
tives of the model equation, i.e., Z — —Xdu/dx = du/dt. This is valid because the 
coefficients satisfy Eqn. (5), which is deduced from the Taylor’s series expansion up 
to the third-order term. On the other hand, the amplification factor of the 3-step 
scheme proposed by Jameson et al., Eqn. (6), can be derived as 

g = l + Z+ l -Z 2 + (24) 

It is obvious that the scheme is not third-order accurate. 

A similar analysis can be applied to the 4-step methods. Identical forms of the 
amplification factors of both 4-step methods of concern (Eqns. (10) and (11)) are 
obtained, such as 


9 = l + Z + lz 2 +l2 3 + ^Z 4 . (25) 

Unlike the case of the 3-step schemes, the effect of the order of accuracy of these 
two 4-step schemes does not appear in the expression of the amplification factor. 
This is because the amplification factor is derived based on the linearized equation. 
Only by using Eqns. (5) and (9), which take into account nonlinear terms, can one 
justify the order of the accuracy of the RK schemes. 

The remaining task is to derive the explicit form of the spatial discretization 
operation, Z, of compact difference schemes. The fourth-order compact difference 
method, Eqn. (12), can be cast into the operator-type by defining 


S U{ — U-l- |_i 2 Ui -|- Ui — (26) 

where Ui could be any flow property of interest at grid point i. As a result, the 
fourth-order method can be rewritten as, 


P\ fdu\ = _L_ 

+ 6 J \dx J • 2 Ax 


(u i+ i - Ui- 1 ) + Q( Ax 4 ). 


(27) 


This equation allows us to express ( du/dx)i in an explicit form, 




To proceed, we substitute this explicit, discretized form, Eqn. (32), into the model 
equation, and we obtain 


Z (4) = 


6Fsin(k)i 
4 + 2 cos(k ) ’ 


(29) 
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where F is the CFL number which is defined as F = \At/Ax , and k is the normal- 
ized wave number (k — 2 irk/ K) . 

It is interesting to note that if the solution reaches a steady state solution, i.e., the 
time derivative term is zero, the operator (1 + 6 2 /6)~ l in the discretized equation 
becomes futile and the spatial discretization is represented by («;+i — ■u;_i)/2A:r, 
which is only second-order accurate. However, the steady state solution of the one- 
dimensional wave equation is a constant and the spatial accuracy is meaningless. 
On the other hand, the accuracy of multi-dimensional calculations is more complex. 
For example, consider a two-dimensional version of the model equation discretized 
by the fourth-order compact difference method, and we have 


du 

m 


+ X a 


1 + -# 

6 


-l 


Uj+lJ Uj — i , j 

2Ax 



Uuj+l 'U’ij — l 

2Ay 


= 0. 


(30) 

Again, the operators (1 +61/6) and (1 + 6^/6) are scalar tridiagonal matrices with 
the dimensions IL x IL and JL x JL, respectively. IL and JL are the numbers 
of the grid nodes in the x and y directions of the computational domain. When 
IL = JL , two operators are identical. We then multiply Eqn. (30) by the operator 
(1 + 61/ 6) and obtain the steady state equation as 


l,jf l,j 


+ A, 




= 0. 


(31) 


2Az y 2A y 

Therefore, the steady state solution is only second-order accurate. 

Similarly, the spatial discretization of the sixth-order compact difference scheme 
can be represented in the operator-type such as, 


du 

dx 



-l 


1 

60Ax 


2 T 28«i + i 


26m-i ~ Ui- 2 ) + 0( Ax 6 ). 


(32) 


And we obtain, 


Z (6) = - 


F[A sin(A;) cos (k) + 56 sin(A:)]z 
12[2 cos (k) + 3] 


(33) 


Again, when solving the steady state solution with same number of grid nodes in 
each direction of the computational domain, the spatial discretization is represented 
by (ui + 2+28ui + i~28ui-\ —Ui-2,)/66Ax. Unlike the case of the fourth-order scheme, 
this representation does not match any conventional central difference scheme and 
it is at most second-order accurate. 


According to the above discussion, we can obtain the amplification factor g(k) for 
various combinations of the Runge-Kutta methods and compact difference schemes. 
The amplification factor g(k) is a complex number and can be expressed as g(k) — 
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exp[iCo(k)], where Co is the normalized frequency and is defined as w = 27rtuAf/T, 
uj is the frequency, r is the the time period of function uf, and the phase speed 
(A) is equal to L/r. Here, the dissipative and dispersive artifacts of the numerical 
schemes can be assessed: 

1. Dissipation. The normalized frequency u>is a complex number {Co = a + i/3) and 
its imaginary part represents the magnitude of the amplification factor, i.e., 

g(k) = = e"<V“ 

wbi = e ~ 0 

The magnitude of the amplification factor is the artificial dissipation. When 
\9\ > 1, the scheme is unstable. For the calculations of unsteady flows, we want 
|g| to be less than unity but very close to it to ensure numerical stability with 
minimum artificial dissipation. In the following section, we plot |<?| against k. to 
illustrate the artificial dissipation. 

2. Dispersion. According to Eqn. (34), the relation of a(k) and k represents the 
artificial dispersion. We plot a = a/F against k to show phase velocities. Notice 
that the model equation is dispersionless and the phase velocity is a constant, i.e., 
A. After being normalized by the CFL number, the exact solution is a straight 
line with 45° angle on the plot of a against k. 

Figure 1 shows the results of the Fourier analysis of the third-order Runge-Kutta 
(RK3) method combined with fourth (CD4), and sixth (CD6) order compact differ- 
ence schemes and the conventional second-order central difference scheme (CD2). 
The figures show the dissipative as well as dispersive effects at CFL numbers from 
0.4 to 1.4 with an increment of 0.2 between neighboring curves. Figure la shows 
the dissipation of the RK3-CD6 scheme. The method is unstable for CFL numbers 
greater than 0.8. As the order of spatial differencing decreases (compare Figs, la, 
lc, and le), the limit of the CFL number increases for stable calculation. 

Figure lb shows the dispersive effect of the RK3-CD6 scheme. For a CFL number 
of 0.4, the phase velocity is correct for wave numbers up to 1.8. The phase velocities 
are slower than they should be at large wave numbers. Increasing the CFL number 
makes phase velocities deviate from the 45° straight line at smaller wave numbers. 
Decreasing the CFL number merges the curves together to reach an asymptotic 
curve. However, the dispersive effect at high wave numbers does not improve. 
Comparing Figs, lb, Id, and If shows that the increase of the order of the spatial 
differencing reduces the numerical dispersion at high wave numbers. Specifically, a 
significant improvement is achieved by changing the spatial differencing from CD2 
to CD4, whereas only a limited gain is obtained by switching from CD4 to CD6. 

Figure 2 shows the results of the Fourier analyses of the fourth-order Runge-Kutta 
(RK4) method combined with various central difference schemes. The CFL numbers 
are the same as that in Fig. 1. Similar to the case of RK3, for the same CFL number 
and wave number, e.g., CFL = 1.4, Co = 1.5, higher-order spatial discretization 
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introduces more artificial damping (see Figs. 2a, 2c, and 2e) and therefore reduces 
the CFL number limit for stable calculation. Again, the dispersive error at high 
wave numbers decreases as the order of the spatial differencing increases (see Figs. 
2b, 2d, and 2f). 

Figures, lc and 2c can be compared to show the difference of the dissipation 
effects between the RK3 and RK4 methods. For the same CFL and wave numbers, 
the RK4 method introduces more artificial damping, and a larger CFL number 
could be used. On the other hand, Figs. Id and 2d show that an increase of the 
order of the time marching scheme does not improve the dispersive effect at high 
wave numbers. 

From the above discussion, it is clear that reliable solutions of the finite difference 
schemes are at low wave numbers. For example, for the RK4-CD6 method at CFL 
= 0.8 (see Figs. 2a and 2b) the solution with wave numbers less than 1/37T (6 grid 
nodes per wave) is fairly accurate. Numerical solutions with higher wave numbers 
(wave length less than 6 grid nodes) suffer significant dispersive and dissipative 
errors. On the other hand, for the conventional RK4-CD2 method at the same 
CFL, 12 to 16 grid nodes per wave are needed for an accurate solution. 

It is interesting to note that compact difference schemes have no dissipative effect 
at the highest wave number resolved by a given numerical grid, i.e., two grid nodes 
per wave (see Figs, la, lc, le, 2a, 2c, and 2e). Nevertheless, a significant dispersive 
error is introduced to these highest-wave-number waves and cause the even-odd de- 
coupling of the numerical solutions. Furthermore, applying the compact difference 
scheme twice to calculate the viscous terms of the Navier Stokes equations does not 
eliminate the erroneous oscillation, owing to the linearity of the operation. These 
high-wave-number waves continue oscillating with erroneous phase speeds through- 
out the course of computation and eventually destroy the solution. It is therefore 
appropriate to impose a small amount of high-order artificial damping to filter out 
these waves while at the same time keeping the resolution at low wave modes intact. 
Figure 3 shows the dissipation and dispersion effects of the RK4-CD6 method at 
CFL=0.8 with various amounts of sixth-order artificial damping, defined as 

JJ 

A.D. = — + U-i — 3 — 6(uj+ 2 + Ui— 2 ) -(- 15(w;+3 + 3) — 20 Uj], (35) 

8 

The range of 77 is from 0.01 to 0.05 with an increment of 0.01 between the neigh- 
boring curves. Comparison between Fig. 3 and Figs. 2a and 2b shows that no 
additional damping at low wave numbers is introduced into the system by the im- 
posed artificial dissipation for 77 < 0.03, whereas the undesirable high-wave-number 
waves are dissipated. 

2.4 Numerical Examples 

2.4.1 Acoustic Admittance of A Nozzle Flow 

The first case is a forced oscillatory quasi-one-dimensional flow in a converging 
nozzle. The governing equations are 
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where 


dQ. dE 

dt dx 


H 


Q=(/^a, E=(pJ + p\a, H=(pgV 

\ e J \(e + p)uj V 0 / 


(36) 


(37) 


p is density, p is pressure, and e is the total energy defined as e = p(C v T +\u 2 ). C v 
is the constant volume specific heat. The variable a is the cross sectional area and 
is prescribed as a function of x. The theoretical solution of the acoustic admittance 
of a choked nozzle was provided by Tsien 7 under the assumption that the velocity 
of the base flow is a linear function of axial location as shown in Fig. 4. The nozzle 
shape can be inversely derived according to Tsien’s assumption, and we have 






7_1 ,\ *(■»-*> 

7 + 1 / 


(38) 


where 7 is the specific heat ratio and M is the Mach number which can be expressed 
as 



7 — 1 / x \ 2 

2 VW * 


(39) 


The superscript * denotes the property at the nozzle throat. According to Tsien’s 
derivation, the linearized quasi-one-dimensional equations can be manipulated to 
the following form under the isentropic condition, 


2(1 _,)*!£_ 2 ( 1+ 


dz 2 


1+7 / dz 


V- 0 £±Mp = o 


2(7 + 1) 


dP 


(7 + 1)(1 - *)5- - ( 7 - 1 + W)P + (2 + *P)U = 0 
az 


(40) 

(41) 


where 


^ = P(z)e i(3T , 
IV 

- - U(z)e il3T . 
u 


(42) 


u and p are the velocity and pressure of the base flow, /3 is the normalized frequency 
which is defined as fi = u;(l — z)/(c* - u ), and r is the non-dimensionalized time 
which is defined as t — c*t/x*. The independent variable z can be expressed 
in different forms due to the linearity between the base flow velocity u and axial 
location x, and we have 
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z 


X 

X* 

u 

c* 

(7 + 1)M 2 
2 + (7 - 1 )M 2 


(43) 


With Eqn. (43), it is clear that P and U are functions of the Mach number (M) 
with prescribed frequency (/?). Equation (40) is a hypergeometric equation 14 that 
can be solved by a power series expansion. The converged solution does not exist 
in the supersonic region because the Mach number is greater than unity. U(z ) can 
be easily solved with P(z) known as shown in Eqn. (41). Finally, the acoustic 
admittance function defined as A(z) = U(z)/P(z) can be obtained as a function of 
the Mach number. 


In what follows, the procedure of the CFD calculation to compare with Tsien’s 
solution is illustrated. First, the base flow field is obtained by solving the quasi- 
one-dimensional equations, Eqns. (36) and (37), using the RK4-CD2 method with 
the nozzle area ratio prescribed by Eqns. (38) and (39). The results are checked 
by the area Mach number relation 15 and the solution is accurate up to five decimal 
digits. The perturbation at the inlet is obtained by specifying sinusoidal pressure 
fluctuations in terms of magnitude and frequency. With the prescribed pressure 
and isentropic correlation, the temperature fluctuation is also determined. Numer- 
ically, these boundary conditions are enforced by defining a vector k = k(Q) at the 
upstream boundary, such as 


k = 




(44) 


where £1 and £2 are the specified values of p and T. To proceed, Equation (44) is 
linearized to become a function of AQ, such as 

k n+1 =k” + — AQ, (45) 

where k n+1 is equal to the specified pressure and temperature at the time step 
n+1 and dk/dQ is a 3 x 3 matrix. In order to close the system, the null entry 
in the vector k may be filled by the out-running characteristic relation deduced 
from the flow equations. Numerically, the similarity transformation is applied to 
the discretized flow equations (see Eqn. (4)), and we get 


LM - 1 (Q' - Q n ) = 

k= 1 


* = !,•••, N 


(46) 
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where i = 1, • • • , N represents the N-step RK method. Here, M 1 is the eigenvector 
matrix of Jacobian matrix A = cdE/dQ, and L is a selection matrix with zeros and 
ones on the diagonal in such a fashion that the proper outrunning characteristics 
are selected. By combining the imposed conditions, Eqn. (45), with the outrunning 
characteristic relations, Eqn. (46), we form the complete equation at the boundary 
point as, 


(Lr' + |) (Q‘-Q”) =LM-‘Ai£a. l R t - 1 +k” + 1 -k", ^ 

i = ,N 

For the supersonic out-flow condition, Eqn. (46) is used with the selection matrix L 
equal to an identity matrix. In both cases, the out-running characteristic equations 
are solved with one-sided difference as shown in Eqn. (14). In other words, the 
characteristic boundary conditions are always discretized by an upwinding scheme 
which is physically sound and the numerical stability is enhanced. These boundary 
conditions are applied at each of the Runge-Kutta stages. 

The acoustic admittance is a complex number and can be written as A = \A\e' e . 
In the present paper, a small pressure perturbation of 1.1% ( p ' = O.Ollp) is imposed 
at the nozzle inlet. The length of the converging part of the nozzle is 0.9 L* (see 
Fig. 4) and the inlet Mach number is about 0.09. The frequency of the perturbation 
is set at (3 = 6, which corresponds to about 2000 Hz. 

Figure 5 shows the comparisons between the CFD results of the RK4-CD6 method 
and the theoretical solution of the acoustic admittance in terms of the magnitude 
\A\ and the phase angle 6 in the subsonic region of the nozzle. Both the magnitude 
and the phase angle of the acoustic admittance decrease as the flow speeds up. As 
shown in the figure, perfect agreement is obtained for the comparison of |A|, while 
the predicted phase angle is slightly off due to the resolution of the numerical grid 
for the phase angle. In this case, the harmonic content of the solution is limited to 
one frequency with a wave length comparable to the computational domain which 
is resolved by 61 grid nodes. Therefore, all numerical schemes of concern provide 
accurate solutions. The numerical errors of |A| and 6 are tabulated in Table 2. 
There is slight advantage in using the higher-order schemes for the prediction of 
\A\; however, no obvious advantage of using the higher-order scheme for the phase 
angle calculation is observed. 

2.4.2 Shocked Sound Waves 

The second case is the propagation of shocked sound waves in a tube with a 
periodic boundary condition. The governing equations are the same as in the first 
case, namely, Eqns. (36) and (37), with cross section area (a) equal to a constant. 
This case is interesting for its complex harmonic content compared to the first 
case. In addition, the capability of the high-order compact difference schemes for 
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shock capturing can also be studied. At time equal to zero, a sinusoidal pressure 
distribution is given. Because of the periodic boundary condition, only one cycle 
resolved by 61 grid nodes is imposed in the computational domain. According to 
the isentropic condition, the distributions of temperature, density, and speed of 
sound are also determined. The velocity profile is determined by the simple wave 
correlation 16 , such as 


where the average flow properties are denoted by a bar. With the simple wave 
correlation, the wave forms of all flow properties are in phase and the initial con- 
dition of the present CFD computation matches the theoretical analysis provided 
by Morse and Ingard 8 . It is interesting to note that the simple wave correlation is 
an extension of a linear, plane, acoustic wave. For a variation of pressure less than 
5%, the plane wave relations could be adopted, such as 


where p'(x) is the prescribed pressure fluctuation. As shown in Eqns. (48) and (49), 
the wave speeds u + c,u-c and u vary as a result of the flow property distribution. 
The distortion of the wave form is a cumulative effect resulting from the wave speed 
distribution. For simple waves, i.e., all flow properties are in phase, the wave crest 
will quickly overtake the trough and form a shock. 

Figure 6 shows the time history of the pressure fluctuation at one end of the 
computational domain for various finite difference schemes. According to Morse and 
Ingard, the first shock appears after about two cycles for the case of a 10% pressure 
perturbation (p ' /p = 0.1) 8 . All schemes of concern predict the wave steepening rate 
correctly. After the wave shocked, the flow evolution is no longer isentropic and the 
kinetic energy is gradually converted to thermal energy due to the existence of the 
shock wave. As a result, the strength of the shock wave diminishes as time passes. 

The shock front is a combination of many wave modes travelling at the same 
speed. The dispersion error introduced by the finite difference schemes will cause 
the high-wave-number waves to travel with erroneous speeds. As shown in Fig. 6, 
the methods of RK4-CD6 and RK4-CD4 with a small amount of the sixth-order 
artificial damping ( rj = 0.02) crisply resolve the shock except for the over-shoots. 



(48) 



(49) 


, , , sp'(x) 

u{x ) = c{x)— — 
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These over-shoots are caused by the Gibbs phenomenon and can be fixed only by 
TVD type shock-capturing schemes. Almost no difference can be observed between 
the results of the CD4 and CD6 methods. On the other hand, the method of RK4- 
CD6 without background filtering shows that significant high-wave-number waves 
lag behind the shock front because the compact difference scheme introduces no 
dissipative but high dispersive effects on the highest wave number waves. As shown 
in the figure, these oscillations eventually contaminate the whole solution. For the 
conventional RK4-CD2 method, results show significant oscillations of moderate 
wave numbers behind the shock front because of dispersion errors. 

Figure 7 shows the normalized power spectrums of the pressure profiles after 
about 17 cycles calculated by different methods. The analytical solution is plotted 
as the solid line. The power of each wave mode is roughly inversely proportional 
to the square of the wave number (a 1 /n 1 2 ). Since 61 grid nodes are used, only 
30 Fourier modes are resolved for the power spectrum (the other 30 modes are the 
complex conjugates). Clearly, the method of RK4-CD2 has significant errors in low 
wave modes. On the other hand, the methods of RK4-CD4 and RK4-CD6 compare 
well with the analytical solution. 

2.4.3 Vortex Propagation in an Uniform Flow 

A Lamb vortex propagated in an uniform flow is chosen as a two-dimensional 
numerical example. The vortex can be characterized by the circulation T and the 
core radius a. The azimuthal velocity uq at a distance r from the vortex center is 
given as, 


T r 

Ue ~ 2tt r 2 + a? ’ 


(50) 


The flow near the vortex center is a rigid-body rotation (uo a r). The flow far 
outside the core is irrotational {u e <x 1/r) with u e decreasing as r increases. Eqn. 
(50) is a continuous function to connect the two extremes. With the prescribed 
velocity field, the pressure and density distributions of the vortex can be determined 
by the momentum and the energy equations, 


dp _ «f 

dr P r ’ 


1 ' o IL ° ’ 

7 - 1 P 2 


(51) 

(52) 


where h a is the total enthalpy and is set to be a constant such a s h Q = jp/i'r - 1 )p 
with the free stream condition denoted by a bar. To proceed, substitute Eqns. (50) 
and (52) into Eqn. (51) and integrate the equation over r. As a result, the pressure 
distribution is obtained. Consequently, the density distribution and the whole flow 
field is determined. The solutions of this stationary vortex can be superimposed 
to any uniform flow with arbitrary speed. Physically, this process may be inter- 
preted as a stationary vortex being observed from a moving coordinate system with 
constant velocity. Thus, the vortex in an uniform flow can be constructed as, 
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V = V + V 1 , 

where the velocities of the background flow are denoted by a bar and the superscript 
' denotes the vortex velocities specified by Eqn. (50). The pressure and density 
distribution of the moving vortex is the same as that of the stationary vortex and 
may be obtained from the solutions of Eqns. (51) and (52). 

The boundary condition of the present case is an extension of the characteristic 
type treatment discussed in Case 1. Essentially, only one-dimensional character- 
istics (derived from two-dimensional flow equations) normal to the computational 
boundary are considered. For the purposes of this discussion, the subsonic out-flow 
condition is considered. The coupled equations of three out-running characteris- 
tics and one specified boundary condition, similar to that in Eqn. (47), should be 
solved. For steady state calculations, a back pressure ( pi , ) is specified to regulate 
the flow rate, such as 


/°\ 

k = ° (54) 

\P 6/ 

Similarly, the dimension of vectors Q and R is four and the matrix M -1 is a 4 x 4 
eigenvector matrix for the flux vector normal to the computational boundary. 

For a non-reflective boundary condition, Giles’ formulation 17 instead of the back 
pressure is used to fill the entry for the specified boundary condition, such as 


dci 

~dt 


+ (0,u,0,v) 


d_ 

dy 



= 0 


(55) 


where y is in the direction parallel to the computational boundary. The variables 
d,i = 1, • • • , 4 are the characteristic variables and can be obtained by the similarity 
transformation from the non-conservative form equations as illustrated by Giles. In 
our case, the characteristic variables are derived from the conservative-form equa- 
tions using the same eigenvector matrix M -1 as in the aforementioned discussion. 
Giles’ non-reflective formulation, Eqn. (55), is relatively simple to used with an 
existing one-dimensional characteristic boundary condition. Nonetheless, according 
to Giles’ analysis, some two-dimensional effect is considered in the equation. Nu- 
merically, Eqn. (55) may be discretized according to the finite difference scheme 
of the interior nodes and combined with the discretized out-running characteristic 
equations (the two-dimensional version of Eqn. (46)) to form the complete subsonic 
out-flow boundary condition. It is interesting to note, however, according to Huff’s 
study 18 and our experience, that stretching the numerical grid nodes downstream 
to dampen the outgoing unsteady waves is just as effective. 
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For the in-flow conditions, the characteristic-type treatment combined with Giles’ 
equation (different from Eqn. (55)) may be adopted. For the present calculations, 
however, the upstream condition is relatively insensitive to various forms of non- 
reflective treatment as long as the proper out-running characteristic equation is 
selected and solved with the prescribed incoming conditions similar to that in Eqn. 
(47). In the present case, constant total pressure and total temperature are pre- 
scribed as the forcing boundary conditions upstream. 

As in the one-dimensional cases, dissipative and dispersive effects of various 
schemes are assessed. The prescribed vortex flow field contains a broad band of 
frequencies due to the distribution of the azimuthal velocity. Theoretically, all wave 
modes travels at the same speed to ensure the integrity of the vortex structure. 
For numerical methods with dispersive error, the shape of the vortex could deform, 
even break up in the later stage of the time marching procedure. In addition, the 
dissipation effect of finite difference schemes can be evaluated by the conservation 
of the sharp pressure dip at the center of the vortex propagated in the numerical 
grid. 

In the present calculations, the Mach number of the background flow is 0.4. The 
grid size is 301 x 91 in the streamwise and transverse directions, respectively. Uni- 
form grids are used in the axial direction and the transverse grids are stretched near 
the outside boundary. The CFL number calculated based on the background flow 
is about 0.7 for all calculations. As discussed before, a small amount of background 
filtering (77 = 0.02) is applied for all calculations. The core radius (a) is about 1 cm 
and is resolved by about 4 grid nodes. 

Figure 8 shows the vorticity and Mach number contours of the initial condition 
prescribed by Eqns. (50) - (53). Figure 9 shows the contours after the eddy prop- 
agates about 60 core radii downstream as simulated by various numerical schemes. 
The comparison between Fig. 8 and Fig. 9 shows that the structure of the eddy is 
retained by the compact difference schemes (CD4 and CD6). In contrast, the eddy 
predicted by the second-order central difference is shattered due to the excessive 
dispersive error. 

Figure 10 shows pressure distributions of the eddy at various instances. In this 
figure, the x axis represents the streamwise locations non-dimensionalized by the 
core radius of the vortex and the y axis is the pressure. For both CD4 and CD6 
methods, the pressure at the vortex center increases about 1 % through the process. 
In comparison, the results of the second-order scheme show pressure fluctuations 
with an overall increase about 3 %. The pressure fluctuation predicted by the 
second-order scheme is due to the deviation of the vortex path. 

2.4.4 Vortex Pairing 

Finally, the calculation of the single vortex is extended to the simulation of the 
vortex pairing. The vortex pairing is the controlling mechanism for the growth of 
the mixing layer 10 . In theory, vortex pairing occurs when the distance between 
two vortices is less than a threshold value. Unfortunately, no theoretical analysis 
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is available for compressible flows. In the present paper, the RK4-CD6 method is 
used to simulate the pairing process to demonstrate the resolution of the high-order 
compact difference scheme. 

The initial condition is specified by two identical vortices placed 5 core radii apart 
in a quiescent gas. The core radius is 1 cm and circulation is 15 m 2 /s. At the center 
of each vortex, there is a pressure deficit about 15% compared to the ambient gas. 
The grid size is 201 x 201. Uniform grids are used at the center of the computational 
domain to resolve the vortices. The grids are slightly stretched in all four directions 
to prevent erroneous wave reflection. In addition, one-dimensional characteristic 
equations combined with Giles’ unsteady, subsonic, out-flow equation is solved on 
the boundary as the non-reflective boundary condition. 

Figure 11 shows the contours of the vorticity magnitude at various stages of the 
vortex interaction. The whole sequence is about one and a half revolutions. Finally, 
a single larger vortex emerges as the result of the vortex pairing interaction. Figure 
12 shows the corresponding Mach number contours for the same flow. 

3. Concluding Remarks 

In this work, the quasi-one-dimensional and two-dimensional Euler solvers us- 
ing various combinations of the Runge Kutta methods and the compact difference 
schemes were developed for numerical simulations of unsteady flows. The accuracy 
of the finite difference schemes is assessed by Fourier analysis and numerical exam- 
ples in terms of numerical dissipation and dispersion. The dispersive characteristic 
is improved by high-order compact difference schemes compared to the second-order 
central difference. The increase of the order of time stepping scheme, on the other 
hand, enlarges the CFL limit for stable computations. In particular, significant im- 
provement of the dispersive effect is obtained by adopting the fourth-order compact 
scheme (6 to 8 grid nodes to resolve a wave) instead of the conventional second-order 
central difference (12 to 16 grid nodes for one wave). The use of the sixth-order 
compact scheme (5 to 8 grid nodes for one wave), however, gains little improvement 
compared with the fourth-order scheme. It was also found that the compact dif- 
ference schemes have no dissipative but high dispersive effects to the highest-wave- 
number waves resolved by a given numerical grid. Consequently, a small amount of 
background filtering is necessary to dissipate the high-wave-number waves and, at 
the same time, keep the low-wave-number solution intact. Other issues such as the 
order of accuracy of the Runge-Kutta schemes for nonlinear equations are analyzed. 
Specifically, the criteria for the 3 and 4-step methods to be third and fourth-order 
accurate are derived. The accuracy of the compact difference schemes for the steady 
state solution is also addressed. 

In general, simulation of unsteady flow provides an overwhelming amount of in- 
formation. It is our experience that the initial and boundary conditions must be 
carefully set up to obtain interpretable and physically meaningful solutions. For 
practical purposes, Giles’ non-reflective equations combined with one-dimensional 
characteristic equations and their implementation to the present numerical scheme 
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were illustrated in detail. In addition, the initial conditions of the simple wave, 
plane acoustic wave, and the Lamb vortex were also provided. Finally, as illus- 
trated in the numerical examples, for flows of simple harmonic content, e.g., one 
frequency in Case 1, the conventional second-order central difference scheme is ad- 
equate provided enough grid nodes are used to resolve the very wave mode. On 
the other hand, for flows of complex harmonic content, the use of the Runge-Kutta 
method combined high-order compact difference schemes shows crisp resolution of 
unsteady flows. 
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Table 1. The Accuracy of the 3-Step Runge-Kutta Methods. 

A t Expansion 3-step R-K methods 

OR R 

1 R (<*31 + <*32 + <*33 )R 

2 [<* 11<*32 + <* 33(<*21 + &22)}RR' 

3 \(R 2 R" + RR' 2 ) \[oc 2 n a i2 + (<* 2i + oc22) 2 a zz \R 2 R" 

+<*1 1 <*22 <*33 RR ' 2 


Table 2. The Relative Error of the Acoustic Admittance Calculation (%). 


Numerical Schemes Error of |A| 
RK4-CD6 0.45 

RK4-CD4 0.52 

RK4-CD2 1.65 


Error of 6 
3.6 
3.3 
4.1 
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a. Dissipation of RK3-CD6 b. Dispersion of RK3-CD6 




c. Dissipation of RK3-CD4 d. Dispersion of RK3-CD4 



e. Dissipation of RK3-CD2 f. Cisaersion of RK3-CD2 


Fig. 1 Dissipation and dispersion characteristics of the RK3 time-stepping combined with 
various spatial discretization schemes. 
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CFL = 0.4 



a. Dissipation of RK4— CD6 



b. Dispersion of RK4— CD6 


0.4 



c. Dissipation of RK4— CD4 


d. Dispersion of RK4— CD4 



e. Dissipation of RK4— CD2 


f. Dispersion of RK4— CD2 


Fig. 2 Dissipation and dispersion characteristics of the RK4 time-stepping combined with 
various spatial discretization schemes. 
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Fig. 3 Dissipation and dispersion characteristics of the RK4-CD6 method with various 
amount of the sixth-order numerical damping. 
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Fig. 4 Configuration of Tsien’s converging-diverging nozzle [7]. 
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Fig. 5 Acoustic admittance calculation using the RK4-CD6 method for the inlet perturbation 
as (3 = 6. p'/p = 0.011. The solid line is Tsien’s theorem and the triangles are the 
CFD results, (a) Magnitude. (b)Phase angle. 
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Fig. 6 Time histories of the pressure fluctuations of the N-wave calculation at one end of the 
periodic domain by various numerical schemes. 
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Fig. 7 Power spectrum of the pressure distribution after about 17 cycles of the N-wave prop- 
agation by various numerical schemes. 
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Fig. 8 Vorticity and Mach number contours of an analytical Lamb vortex. 
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Fig. 9 Vorticity and Mach number contours of the lamb vortex after travelling about 60 core 
diameters predicted by various numerical schemes. 
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Fig. 10 Vortex pressure profiles at the center line at various instances predicted by different 
numerical schemes. 
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Fig. I [ The vorticity magnitude contours for the vortex pairing. 
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Fig- It 


The vorticity magnitude contours for the vortex pairing. 
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The Mach number contours for the vortex pairing. 
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Appendix B 

CMOTT Biweekly Seminars / Technical Meetings 


The purpose of these seminars is to exchange ideas and opinions on the latest 
developments and current state of turbulence and transition research. The speakers 
are invited from within and without of the NASA LeRC, including foreign speakers. 
The seminars were intended not noly to keep the members informed of the latest 
development of local turbulence and transition modeling research but also to increase 
interactions between group members and other researchers at the NASA LeRc. 

The following is the meeting schedule and the abstract of the semainars during 
the reporting period. 
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CENTER FOR MODELING OF 
TURBULENCE AND TRANSITION 

Date: July 3, 1991 

To: CMOTT Members and SVR and IFMD Staff 

From: William W. Liou (6682) 

Subject: CMOTT Biweekly Meeting 



The following is a tentative schedule for the CMOTT biweekly get-together from July 
10, 1991 to August 28, 1991. 

The presentations will be informal and active participation is expected from the at- 
tendants. Soda and snack will be served in the meetings. These meetings complement the 
CMOTT Seminar Series, which are mainly formal presentations. 

We would also appreciate some contributions from you. Subjects related to either the 
theoretical, experimental or computational aspects of turbulence and transition modeling 
are welcomed. Those who are willing to share their experience in these areas can contact 
me or Dr. T.-H. Shih at 6680 for further arrangement. 

The meeting will start at 4:00 p.m. in Room 228, Sverdrup Building. 


July 10, 1991 J. Lepicovsky (61-6753) 

LDV Measurement of Large Structures in a Tone 
Excited Turbulent Jet 

July 24, 1991 C. R. Wang (5865) 

Computations of Turbulence in a Shock/Turbulent 
Boundary Layer Interaction Flow 


August 7, 1991 A. Hsu (61-6648) 

PDF Turbulence Model and Its Applications 


August 28, 1991 C. Steffen (8508) 

DTNS: An Accurate and Efficient Testbed for 
Incompressible Flow Turbulence Modeling 
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Date: 

September 4, 1991 

To: 

CMOTT Members and SVR and IFMD Staff 

From: 

William W. Liou (6682) 

Subject: 

CMOTT Biweekly Meeting 


The following is a tentative schedule for the CMOTT biweekly get-together from 

September 11, 1991 to October 23, 1991 . 

The presentations will be informal and active participation is expected from the at- 
tendants. Soda and snack will be served in the meetings. These meetings complement the 
CMOTT Seminar Series, which are mainly formal presentations. 

We would also appreciate some contributions from you. Subjects related to either the 
theoretical, experimental or computational aspects of turbulence and transition modeling 
axe welcomed. Those who are willing to share their experience in these areas can contact 
me or Dr. T.-H. Shih at 6680 for further arrangement. 

The meeting will start at 4:00 p.m. in Room 228, Sverdrup Building. 

Sept. 11, 1991 T. Bui (5639) 

Implementation of the Chien Low-Re k-e Models into the 
Proteus Code 

Sept. 25, 1991 IC. Ahn (5965) 

A 2-D Oscillating Flow Analysis Using Quasi-steady 
Turbulence Model 

October 9, 1991 J. Schwab (8446) 

Variable-density Turbulence Modeling for Turbomachinery 

October 23, 1991 M. Mawid (5965) 

Multiphase Turbulent Combustion 
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Date: 

November 4, 1991 

To: 

CMOTT Members and SVR and IFMD Staff 

From: 

William W. Liou (3-6682) 

Subject: 

CMOTT Biweekly Meeting 


The following is a tentative schedule for the CMOTT biweekly get-together from 

November 6, 1991 to December 18, 1991. 

This will be the last session of the CMOTT group- meet ings/inforinal- seminars this 
year but the series will resume in mid- January 1992. Thank you for 3'our patience and 
participation through out the year. The group-meetings/ informal-seminars of CMOTT 
are meant not only to serve CMOTT members but also to provide an informal forum for 
those who are involved in transition/turbulence predictions. I thank all the speakers and 
participants who have made these objectives “realizable”. Now, we arc planing for 1992. 
If you have any suggestions or like to give a talk or two in the coming year, please call me 
or Dr. T. H. Shill at 3-6680. In the mean time, don’t forget to imirk your calendar for the 
following talks ! 


HAPPY HOLIDAYS !!! 

The meeting will start at 4:00 p.m. in Room 228, Sverdrup Building. 


Nov. 6, 1991 W. Liou (3-66S2) 

Weakly Nonlinear Models for Turbulent Free Shear Flows (2) 
- A Self-Contained Energy Transfer Model 

Nov 20, 1991 D. Ashpis (3-8317) 

DNS of Disturbances in Boundary Layer Flow 

Dec. 4, 1991 B. Rubinstein (61-6612) 

Analytical Theory of Turbulence and Turbulence Modeling: 
TSDIA and RNG 

Dec. 18, 1991 B. Duncan (61-2998) 

A New Three-Equation Model for Turbulence 
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Date: 


January 30, 1992 


To: 


CMOTT Members and SVR and IFMD Staff 


From: 


William W. Liou (3-6682) 


Subject: CMOTT Biweekly Meeting 

The following is a tentative program for the CMOTT biweekly get-together/seminar from 
February 5, 1992 to March 18, 1992. Ybu are welcomed to join us. 

Thanks to the your suggestions, we have made a few changes from last year’s format. First, 
we have scheduled the CMOTT Seminar Series, which are mainly formal presentations, into 
the biweekly time frame of the CMOTT group-meeting/informal-talks. Also, the abstract 
of each presentation, formal or informal, will be distributed about one week prior to its 
scheduled date. Again, if you are interested in giving a presentation, please contact us. 

The meeting will run from 1:30-2:15 PM in Room 145, Sverdrup Building, unless otherwise 


noted. 


(1) Feb. 5, 1992 D. Davis 


Weakly Nonlinear Vortex/Wave Interactions in 
Incompressible Cross-flow Boundary Layers in Transition 


(2) Feb. 19, 1992 Z. Yang 

A Modeling of Bypass Transition 

(3) March 4, 1992 K. Zaman 


Effect of ’’Delta- Tabs” on the Evolution of Axisymmetric 
Jets 


(4) March 18, 1992 Professor R. M. C. So, Arizona State University 

Near Wall Heat Transfer Modeling 
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Date: March 26, 1992 

To: CMOTT Members and SVR and IFMD Staff 

From: William W. Liou (3-6682) 

Subject: CMOTT Biweekly Meeting 

The following is a tentative program for the CMOTT biweekly get-together/seminar from 
April 1, 1992 to May 13, 1992. You are welcomed to join us. Also, if you are interested 
in giving a presentation, please let us know. 

The meeting will run from 1:30-2:15 PM in Room 228, Sverdrup Building, unless otherwise 
noted. 

(5) April 1, 1992 J. Van der Vegt 

The Development of an ENO-Osher Scheme for Direct 
Simulation of Compressible Flows 

(6) April 15, 1992 J. Goodrich 

Unsteady Time Asymptotic State: Incompressible Results, 

New Directions for Algorithms 

(7) April 29, 1992 T.-H. Shih 

Kolmogorov Behavior of Near- Wall Turbulence and 
Its Application in Turbulence Modeling 

(8) May 13, 1992 Z. Yang 

A Modeling of Bypass Transition 
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Date: June 1, 1992 

To: CMOTT Members and SVR and IFMD Staff 

From: William W. Liou (3-6682) 

Subject: CMOTT Biweekly Meeting 



The following is a tentative program for the CMOTT biweekly get-together/seminar from 
June 10, 1992 to July 22, 1992. You are welcomed to join us. 

The talks will be informal and active participation will be expected from the audience. 

Also, if you are interested in giving a presentation about the progress or some results of 
your own work on turbulence or transition, please let us know. 

The meeting will run from 1:30-2:15 PM in Room 228, Sverdrup Building, unless otherwise 
noted. 


(9) June 10, 1992 J. Zhu 

Finite Volume Computations in Incompressible Flows 
with Complex Geometries 

(10) June 24, 1992 J. Lee 

RPLUS Code and Standard k — e Models and Applications 

(11) July 8, 1992 R. Mankbadi 

Unsteady Turbulent Flows 

(12) July 22, 1992 D. Rigby 

The Effect of Spanwise Variations in Momentum on 
Leading Edge Heat Transfer 
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Biweekly Meeting Series (1) 


Weakly Nonlinear Vortex/Wave Interactions in 
Incompressible Crossflow Boundary Layers in 

Transition 


The instability of an incompressible three-dimensional boundary layer 
is considered theoretically and computationally in the context of vor- 
tex/wave interactions. Specifically the work centres on two low-amplitude, 
lower-branch Tollmien-Schlichting (TS) waves which mutually interact 
to induce a weak longitudinal vortex flow;the vortex motion, in turn, gives 
rise to significant wave-modulation via wall-shear forcing.The character- 
istic Reynolds number is taken as a large parameter and, as a conse- 
quence.the TS waves and the vortex are governed primarily by triple- 
deck theory.The nonlinear interaction is captured by a viscous partial- 
differential system for the vortex coupled with a pair of amplitude equa- 
tions for the wave pressures. Computations were performed for relatively 
small crossflow values.Three distinct possibilities were found to emerge 
for the nonlinear behaviour of the flow solution downstream - an alge- 
braic finite-distance singularity.far- downstream decay or repeated oscil- 
lations - depending on the various parameter values, the input amplitudes 
and the wave angles. 


by 

Dominic Davis 
ICOMP 


Wed., 5 February, 1992 
1:30-2:15 PM 
Room 145, SVR Building 


ABSTRACT 
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A k-e Model for Near Wall Turbulence and 
its Application in Turbulent Boundary Layer 

with/without Pressure Gradient 

by 

Z. Yang 

ICOMP 



Wed., 19 February, 1992 
1:30-2:15 PM 
Room 145, SVR Building 


ABSTRACT 

A k-e model is proposed for turbulent wall bounded flows. In this 
model, turbulent velocity scale and turbulent time scale are used to 
define the eddy viscosity. The time scale used is bounded from below 
by the Kolmogorov time scale. The dissipation equation is reformulated 
using this time scale, removing the singularity of the high Reynolds 
number k-e equation at the wall and rendering the introduction of 
the pseudo-dissipation unnecessary. The damping function used in the 
eddy viscosity is chosen to be a function of R y = - y instead of y + . 
Thus, the model could be used for flows with separation. The model 
constants used are the same as the model constants in the commonly 
used high turbulent Reynolds standard k-e model. Thus, the proposed 
model would reduce to the standard k — e model when it is far away 
from the wall. Boundary layer flows at zero pressure gradient, favorable 
pressure gradient, adverse pressure gradient and increasingly adverse 
pressure gradient are calculated respectively. The comparisons of model 
predictions and the available experimental data are found to be good. 
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Effect of Tabs on the Evolution of Axisymmetric Jets 

by 


Khairul Zaman 

IFMD 


Wed., 4 March, 1992 
1:30-2:15 PM 
Room 145, SVR Building 


ABSTRACT 


Vortex generators, in the form of small tabs at the nozzle exit, can have 
a profound influence on the evolution of an axisymmetric jet. Using 
tabs of certain shapes, the jet cross section can be distorted almost 
arbitrarily. Such distortion is accompanied by elimination of screech 
noise from supersonic jets and a significant increase in jet mixing. Key 
results obtained so far, covering a jet Mach number range of 0.3 and 
1.8, will be summarized in this presentation. Observations will be made 
on the mechanisms of the effect including the likely vorticity dynamics 
in the flow. 
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Near-Wall Modeling of Turbulent Heat Transfer 

by 

Professor Ronald M. C. So 



Mechanical and Aerospace Engineering 
Arizona State University 


Wed., 18 March, 1992 
1:30-2:30 PM 
Room 119, Building 5 


ABSTRACT 


A near-wall two-equation model for turbulent heat fluxes is derived from 
the temperature variance and its dissipation-rate equations and the as- 
sumption of gradient transport. The near-wall asymptotics of each term 
in the exact equations are examined and used to derive near-wall cor- 
rection functions that render the modeled equations consistent with 
these behavior. Thus modeled, the equations are used to calculate 
fully-developed pipe and channel flows with hear transfer. It is found 
that the proposed two-equation model yields asymptotically correct near- 
wall behavior for the normal heat flux, the temperature variance and its 
near-wall budget and correct limiting wall values for these properties 
compared to direct simulation data and measurements obtained under 
different wall boundary conditions. 


CONTACT: T.-H. Shih, PABX 3-5698 
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The Development of an ENO-Osher Scheme for 
Direct Simulation of Compressible Flows 

by 

Jaap Van der Vegt 
ICOMP 


Wed., April 1, 1992 
1:30-2:30 PM 
Room 228, SVR Building 


ABSTRACT 

Direct simulation of turbulence and transition in compressible wall 
bounded flows presents an alternative to investigate important physical 
phenomena which are difficult to measure or study otherwise. It also 
provides data useful for turbulence modeling. A new program is be- 
ing developed which solves the three-dimensional compressible Navier- 
Stokes equations using a higher order, fully implicit and time accurate 
ENO scheme together with Osher flux splitting. In this presentation an 
overview will be given of the numerical scheme and several test cases, 
both for supersonic and subsonic flow, will be presented and further 
improvements will be discussed. 
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Unsteady Time Asymptotic States: Incompressible Results 
and New Directions for Algorithms 

by 

John Goodrich 

IFMD 

Wed., April 15, 1992 
1:30-2:15 PM 
Room 228, SVR Building 

ABSTRACT 

Unsteady time asymptotic flow states for high Reynolds number viscous incom- 
pressible flow problems are presented. Discrete frequency flows are shown for 
the square driven cavity, with periodic cases for Re = 9000 and Re = 9600 , and 
with aperiodic cases for Re = 9700 and Re = 10000 . The algorithm for these cal- 
culations is based on the fourth order PDE for incompressible fluid flow which 
uses the streamfunction as the only dependent variable, it is second order ac- 
curate in space and time, and it has a stability constraint of CFL < 1 . The 
algorithm is extremely robust with respect to Reynolds number. 

The direct numerical simulation of transition and turbulence requires nu- 
merical methods to be more than second order accurate in order to accurately 
represent the relevant scales of the physical processes. Recently developed 
finite difference algorithms are presented for unsteady convection equations, 
including the advection and inviscid Burgers equation in one space dimension, 
and the wave equation treated as a system, with remarks on diffusion equa- 
tions and extension to higher space dimensions. The new algorithms that will 
be discussed all use local stencils, they range from third to seventh order in 
accuracy, they all have the same order of accuracy in both space and time, and 
they are all one step explicit methods (except for diffusion equations). Since all 
of the algorithms use a small local stencil, the number of degrees of freedom of 
known data required for higher order accuracy is obtained by higher information 
density than just the solution data. The use of a two point stencil (for some 
of the methods) allows for arbitrary grid spacing, though a convective stability 
constraint must be observed at each grid point. The use of local data for an 
explicit algorithm with high order accuracy avoids the requirement of using a 
global solution method such as compact differencing or spline based algorithms. 

There will be computational results for all of the algorithms that are presented. 
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Kolmogorov Behavior of Near-Wall Turbulence 
and Its Application in Turbulence Modeling 

by 

Tsan-Hsing Shih 

ICOMP 


Wed., April 29, 1992 
1:30-2:15 PM 
Room 228, SVR Building 


ABSTRACT 

The near-wall behavior of turbulence is re-examined in a way different from 
that proposed by Hanjalic and LaunderM and followersl a M s M<M*l. It is shown 
that at a certain distance from the wall, all energetic large eddies will reduce 
to Kolmogorov eddies (the smallest eddies in turbulence). All the important 
wall parameters, such as friction velocity, viscous length scale, and mean strain 
rate at the wall, are characterized by Kolmogorov microscales. According to this 
Kolmogorov behavior of near-wall turbulence, the turbulence quantities, such 
as turbulent kinetic energy, dissipation rate, etc. at the location where the 
large eddies become “ Kolmogorov " eddies, can be estimated by using both di- 
rect numerical simulation (DNS) data and asymptotic analysis of near-wall 
turbulence. This information will provide useful boundary conditions for the 
turbulent transport equations. As an example, the concept is incorporated in 
the standard k-e model which is then applied to channel and boundary layer 
flows. Using appropriate boundary conditions (based on Kolmogorov behavior 
of near-wall turbulence), there is no need for any wall-modification to the k- 
e equations (including model constants). Results compare very well with the 
DNS and experimental data. 
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A Modeling of Bypass Transition 

by 

Z. Yang 

ICOMP 

Wed., May 13, 1992 
1:30-2:15 PM 
Room 228, SVR Building 


ABSTRACT 

A model for the calculation of bypass transitional boundary layers due to the 
freestream turbulence is proposed. The model combines a near wall k-e model 
proposed for the fully developed turbulent flows with the intermittency of the 
transitional boundary layers. The intermittency factor is assumed to be a func- 
tion of both the free stream turbulence and the shape factor of the boundary 
layer. Transitional boundary layers over a flat plate with different freestream 
turbulence level are calculated using the proposed model. It is found that the 
model calculations agree well with the experimental data and give a better pre- 
diction compared with other low Reynolds number k — e models, which do not 
incorporate the intermittency effect. 
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Finite-Volume Computations of Incompressible 
Flows with Complex Geometries 


by 


J. Zhu 

ICOMP 


Wed., June 10, 1992 
1:30-2:30 PM 
Room 228, SYR Building 


ABSTRACT 

A brief review is given of finite-volume procedures developed at the 
Institute for Hydromechanics, University of Karlsruhe, for calculating 
incompressible elliptic flows with complex boundaries. The procedures 
include: numerical grid generation, higher-order bounded convection 

schemes, zonal solution, simulation of two-phase flows, and near-wall 
turbulence modelling. Various application examples will be given. 
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Biweekly Meeting Series (1992-10) 


Development of the RPLUS code with 
Standard k-e model and Their Applications 

by 

Jinho Lee 

Sverdrup Technology, Inc. 


Wed., June 24, 1992 
1:30-2:30 PM 
Room 228, SVR Building 


ABSTRACT 

The primary goal of this research effort is to develop a CFD tool which can be used in 
a variety of practical supersonic/hypersonic propulsion device development /analysis envi- 
ronments. One focus of this work has been to develop and validate the Reactive Propulsion 
code based on LU Scheme( RPLUS). This effort also includes the development of turbu- 
lence models which can be used in the predictions of highly complex flow environments 
inside of combustors. 

This presentation will cover only a small part of a larger development effort and focus 
will primarily on the analysis, implementation, and development of the turbulence model 
capabilities of the RPLUS code. 

Some of the issues which will be covered are; formulation of the turbulence models, 
the numerical technique used to solve the turbulence model equations, and modeling of 
compressibility effects. The primary numerical technique used in the RPLUS code is 
the LU-SSOR(LU scheme based on Successive Symmetric Over Relaxation) technique. 
Therefore, the turbulent kinetic energy transport and dissipation transport equations are 
also solved using the LU-SSOR numerical technique. 

Both two and three dimensional turbulence model development are being developed 
for the RPLUS code. However, the majority of the presentation will focus on the devel- 
opment of the two dimensional k-e models for the RPLUS code. Issues regarding com- 
pressible wall-function boundary conditions and compressibility effects will be addressed. 
Both low and high Reynolds number forms of the k-e models are being developed. The 
“standard low Reynolds number model of Launder- Sharma <md Chien has been used in 
this study. The problems of primary interests are supersonic turbulent boundary-layers, 
shock-wave/bounda^-laj'er interactions, and shear-layers in two or three dimensional en- 
vironments. 
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Unsteady Turbulent Flows 

by 

Reda R. Mankbadi 

NASA Lewis Research Center 


Wed., July 8, 1992 
1:30-2:30 PM 
Room 228, SVR Building 

ABSTRACT 

Current research activities emphasize computation/modelling of turbu- 
lent flows when basic flow is time-periodic. Wall-bounded flows and free 
shear flows exhibit different features when the basic flow is unsteady; 
and as such, different approaches are used to model them. 

(A) In wall-bounded, oscillatory flows, two approaches are used to 
model turbulence: (I) Turbulence is assumed to behave in a quasi-steady 
manner, and steady-state models are directly extended to the unsteady 
case. This approach fails at high frequencies of oscillations. (II) Rapid 
distortion theory (RDT) is successfully adapted to aid in turbulence 
modelling of highly unsteady flows (high frequencies). The eddy viscosity 
hypothesis is replaced by the ratio of turbulent stresses/kinetic energy; 
which is given by RDT as a function of the accumulated rate of strain. 

(B) In free shear flows (naturally unsteady, or excited to be un- 
steady), two approaches are investigated: (I) The large-scale (organized, 
coherent) component is modelled as instability waves interacting with 
each other as well as with the mean flow and the fine-scale (random, 
background) turbulence. Integrated kinetic-energy equations are then 
obtained for each scale of motion. The approach is successful in pre- 
dicting results in good agreements with experiments in which excitation 
devices are used to control jet mixing and turbulence. (II) The other 
approach adopted is Large-Eddy Simulations (LES) with application to 
predicting the far-field noise of a supersonic jet. 

CONTACT: William W. Liou, PABX 3-GG82 
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The Effect of Spanwise Variations in 
Momentum on Leading Edge Heat Transfer 


A study of the effect of spanwise variation in momentum on leading 
edge heat transfer is discussed. Numerical and experimental results 
are presented for a circular leading edge and for a 3:1 elliptical leading 
edge. Direct comparison of the two-dimensional results, that is with no 
spanwise variations, to the analytical results of Frossling is very good. 
The numerical calculation, using the PARC3D code, solves the three- 
dimensional Navier-Stokes equations, assuming steady laminar flow on 
the leading edge region. Experimentally, increases in spanwise averaged 
heat transfer coefficient as high as 50% above the two-dimensional value 
were observed. Numerically, the heat transfer coefficient was seen to 
increase by as much as 25% percent. In general, the circular leading 
edge, under the same flow conditions, produced a higher heat transfer 
rate than the elliptical leading edge. As a percentage of the respective 
two-dimensional values, the circular and elliptical leading edges showed 
similar sensitivity to spanwise variations in momentum. By equating the 
root mean square of the amplitude of the spanwise variation in momen- 
tum to the turbulence intensity, a qualitative comparison between the 
present work and turbulent results was possible. 


by 

David Rigby 

Sverdrup Tech. INc. 


Wed., July 22, 1992 
1:30-2:30 PM 
Room 228, SVR Building 


ABSTRACT 


CONTACT: William W. Liou, PABX 3-6682 
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